View of LowerMatrix is not working with ldiv! and mul!

Specific Domains Numerics
1

Hello,

I am computed a large Cholesky factor L \in R^{N,N}, and I would like to multiply and divide an array of dimension n, M with n \leq N. I would like to use view(L,1:n,1:n) but ldiv!, mul! are not defined for views.

using LinearAlgebra
N = 50
M = 500

X = randn(N, M)

L = cholesky(X*X' + 1.0* I).L

n= 10

Y = randn(n, M)

ldiv!(view(L,1:n,1:n), Y) 

mul!(view(L,1:n,1:n), Y)

The error message for mul!:

 mul!(::AbstractArray, ::AbstractArray, !Matched::Number, !Matched::Number, !Matched::Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/generic.jl:126
  mul!(::AbstractArray{T,2} where T, ::Union{DenseArray{T,2}, Base.ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T, !Matched::Transpose{#s662,#s661} where #s661<:Diagonal where #s662, !Matched::Number, !Matched::Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:341
  mul!(::AbstractArray{T,2} where T, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, !Matched::Transpose{#s662,#s661} where #s661<:(Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) where #s662, !Matched::Number, !Matched::Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:345
2

This problem is what ArrayLayouts.jl is designed to solve, but note that view(L,1:n,1:n) does not know from the type information that the view is lower triangular: this means its not clear that ldiv! can be overwrite the result. view(L,Base.OneTo(n),Base.OneTo(n)) is a bit better though not yet supported, but would be easy to add a LowerTrapezoidalLayout for this case so that ArrayLayouts.ldiv!(view(L,Base.OneTo(n),Base.OneTo(n)), Y) works.

I don’t have time to do this myself but happy to walk you through it if you are willing to make a PR.

1 Like
3

Thank you @dlfivefifty. Yes, I will be happy to work on LowerTrapezoidalLayout. How can I start?

4

Look at TriangularLayout in ArrayLayouts.jl/src/memorylayouts.jl and triangular.jl

Probably should have said TrapezoidalLayout as the upper/lower is determined by a character.

5

(Please file an issue on github to discuss further questions)

6

As a short term work around, you can reverse the view and the LowerTriangular wrapping, i.e. something like

julia> L = LowerTriangular(randn(4,4));

julia> ldiv!(LowerTriangular(parent(L)[1:3, 1:3]), ones(3))
3-element Array{Float64,1}:
 -0.748637727829292
  4.553737980760209
 -8.60259724047278
1 Like