Can Not Use Triangular Matrix from Cholesky in Any Way

General Usage
1

This is a regression in Julia 1.0. Consider this code

using LinearAlgebra
Umat = [ 4.10229+0.0im -0.00360251-2.02412im ; -0.00360251+2.02412im 1.30966+0.0im]
R = cholesky(Umat) ;
display(R)

V = Array(R) ;
display(V)

X = zeros(Complex{Float64}, 2, 2, 2) ;
X[2,:,:] = R ;

The output of this in Julia 1.0 is:

Cholesky{Complex{Float64},Array{Complex{Float64},2}}
U factor:
2×2 UpperTriangular{Complex{Float64},Array{Complex{Float64},2}}:
2.02541+0.0im -0.00177866-0.999363im
⋅ 0.557612+0.0im
2×2 Array{Complex{Float64},2}:
4.10229+0.0im -0.00360251-2.02412im
-0.00360251+2.02412im 1.30966+0.0im
ERROR: LoadError: MethodError: no method matching setindex_shape_check(::Cholesky{Complex{Float64},Array{Complex{Float64},2}}, ::Int64, ::Int64, ::Int64)
Closest candidates are:
setindex_shape_check(::AbstractArray, ::Integer…) at indices.jl:154
setindex_shape_check(::AbstractArray{#s57,1} where #s57, ::Integer, ::Integer) at indices.jl:196
setindex_shape_check(::AbstractArray{#s57,2} where #s57, ::Integer, ::Integer) at indices.jl:200
…
This code works perfectly in Julia 0.6.4 once cholesky is replaced with chol() and the using LinearAlgebra statement is removed.

So it appears that the Array() constructor fails to work with a triangular matrix as an argument. (It worked before). Furthermore there no longer appears to be a working auto-conversion to a dense matrix when I attempt to load the triangular matrix into a dense array.

Even worse, this fails
convert(Array{ComplexF64,2}, R)
with:

ERROR: MethodError: no method matching Array{Complex{Float64},2}(::Cholesky{Complex{Float64},Array{Complex{Float64},2}})
How am I supposed to convert these matrices to regular dense arrays when I need to?

Furthermore no matrix operations work on triangular matrices namely,
R’ * R
R+R
R * 2

etc.
I believe the fact that convert(Array, R) does not work properly, is definitely a bug, but my issue was immediately closed.

This is a big step back in usability, unless there is yet another way to convert to a dense matrix that is somehow hidden in the documentation.

FYI it appears that the conversion routines suggested revert R back to R’ * R
thus
Matrix ( R ) = R’ * R
How is this supposed to be helpful? There is no way to use the Cholesky factor at all with this convention.

2

The usual response for this type of questions is to run the old code in Julia 0.7 to see the deprecations and try to correct them. If you have a code that works in 0.6.4, and doesn’t run at all in 0.7 without any deprecation warning, then there’s an issue. Most of these issues are due to some big changes in the names of functions and other things like that, I don’t think the functionality disappeared.

Also, try to enclose your code with three backticks (```), that helps readability.

1 Like
3

Rather than

try

R = cholesky(Umat).U

The cholesky function returns a factorisation object from which you can get either the upper triangular part or lower triangular part as desired. See the docs at https://docs.julialang.org/en/latest/stdlib/LinearAlgebra/#LinearAlgebra.cholesky

2 Likes
4

I ran it in 0.7. You get the same problem. All conversions to dense matrices cause upper triangular R to be converted to R’ * R. That has to be a bug.

The loading into the dense subarray has other issues. It claims that it’s broadcasting but it’s not. In fact X[2,:,:] = V does not complain at all, since V is an ordinary dense array.

5

Why is the default conversion back to R’ * R then?
convert(Array, R) = R’ * R essentially.

6

cholesky(A) does not return the upper (or lower) triangular factor. It returns a factorization object, which is just another representation of the original matrix A, so it is only natural that convert(Array, cholesky(A)) returns something that is == to the original matrix.
If you want the upper or lower triangular factor (i.e. what chol did in Julia v0.6) then you should use cholesky(A).U or cholesky(A).L respectively.

2 Likes
7

I’m not sure what you mean. You can use the upper triangular form in most matrix operations directly or convert it into a regular matrix as you want. For example,

julia> Umat = [ 4.10229+0.0im -0.00360251-2.02412im ; -0.00360251+2.02412im 1.30966+0.0im]
2×2 Array{Complex{Float64},2}:
     4.10229+0.0im      -0.00360251-2.02412im
 -0.00360251+2.02412im      1.30966+0.0im    

julia> R = cholesky(Umat)
Cholesky{Complex{Float64},Array{Complex{Float64},2}}
U factor:
2×2 UpperTriangular{Complex{Float64},Array{Complex{Float64},2}}:
 2.02541+0.0im  -0.00177866-0.999363im
         ⋅         0.557612+0.0im     

julia> RU = Array(R.U)
2×2 Array{Complex{Float64},2}:
 2.02541+0.0im  -0.00177866-0.999363im
     0.0+0.0im     0.557612+0.0im     

julia> RU'*RU ≈ Umat
true
1 Like
8

Well you’ve helped me inasmuch as I didn’t notice that the Cholesky factorization also became an object similar to the changes in qr(), thus I have to use the .U operation to get what I want.

However the object itself can be converted to a dense matrix. ie in your example one can do

R = cholesky(Umat)
RU = Array(R)

However RU will be R.U’ * R.U not R. This propagated all sorts of interesting hard to find bugs in my code.

9

R is basically R.U' * R.U and hence also equal to RU. That is what I meant with

so the object is a representation of A, not the upper or lower factor. For example, it would be weird if

cholesky(A)\x

was not equal to

A\x

so in this example it is pretty clear that cholesky(A) is the same as A, but with different underlying storage.

2 Likes
10

I won’t argue that point due to it’s somewhat pedantic nature, but it was certainly unexpected for me. I think it might have been better for the convert to call an error and perhaps point out the existence of the subfields .U and .L. That would have saved me some time.

Also this was not an error that running under 0.7 would have caught as well. I’ll guess I’ll have to be more careful next time. Thanks for the help.

11

Well, the following message is printed:

julia> using LinearAlgebra

julia> A = rand(2,2); A = A'A;

julia> R = chol(A)
┌ Warning: `chol(A::AbstractMatrix)` is deprecated, use `(cholesky(A)).U` instead.
│   caller = top-level scope at none:0
└ @ Core none:0
2×2 UpperTriangular{Float64,Array{Float64,2}}:
 0.829981  0.527958
  ⋅        0.781579
1 Like
12

Well I outsmarted myself because I already knew chol() was deprecated and never saw that LOL.

1 Like
13

Perhaps the manual could mention that all factorizations work this way. Can you make a pull request?

1 Like
14

That would certainly help. However I didn’t find chol() in the v 0.7 deprecated or breaking sections. I looked for it there. It might be a good idea to update the manual there, since that’s where a lot of the surprising changes are first documented. The fact that cholfact becomes cholesky is there, but not mention of chol() itself. I should have read the new cholesky entry more carefully, but I was in a hurry and the wrong conversion made some of my errors go away.

15

It might be a good idea to update the manual there

Submitting a documentation fix is easy because you can do it in your browser: Editing files - GitHub Docs. Give it a shot, @mcbro!

4 Likes