LinearAlgebra fails on Symmetric complex matrices

[Timeroot]

I have a collection of positive-definite, symmetric, complex matrices I need to quickly solve linear systems on. Given the positive definiteness, I initially tried cholesky but that gave me the error:

MethodError: no method matching cholesky(::Symmetric{ComplexF64, Matrix{ComplexF64}})

so I figured that wasn't supported. I tried factorize figuring that would pick the right method for the job. It fell back to Bunch-Kaufman, but then gave me the error:

ArgumentError: adjoint not implemented for complex symmetric matrices

Stacktrace:
 [1] adjoint(B::BunchKaufman{ComplexF64, Matrix{ComplexF64}})
   @ LinearAlgebra /opt/julia-1.7.1/share/julia/stdlib/v1.7/LinearAlgebra/src/bunchkaufman.jl:284

Ultimately the only current workaround is to explicitly turn my Symmetric{ComplexF64} into a Matrix{ComplexF64} and then call factorize, or to directly call lu myself. Either way, this loses any optimizations that could be done using the symmetry or, even better, the positive definiteness.

Besides that BunchKaufman should work for complex symmetric matrices, Cholesky decomposition also works for PSD complex symmetric matrices; the difference is that M == L * transpose(L) instead of L * L'. Here's an implementation I made of a modified Cholesky decomposition for complex symmetric matrices:

function my_chol!(A::AbstractMatrix, ::Type{LowerTriangular})
    LinearAlgebra.require_one_based_indexing(A)
    n = LinearAlgebra.checksquare(A)
    @inbounds begin
        for k = 1:n
            Akk = A[k,k]
            for i = 1:k - 1
                Akk -= A[k,i]*A[k,i]
            end
            A[k,k] = Akk
            Akk, info = sqrt(Akk), 0
            if info != 0
                return LowerTriangular(A), info
            end
            A[k,k] = Akk
            AkkInv = inv(Akk)
            for j = 1:k - 1
                @simd for i = k + 1:n
                    A[i,k] -= A[i,j]*A[k,j]
                end
            end
            for i = k + 1:n
                A[i,k] *= AkkInv
            end
        end
     end
    return LowerTriangular(A), 0
end

(obviously the above can be cleaned up, I left it as-is to maximize similarity to existing chol! code.)

If you want to try this out, here's a snippet with two simple test cases to try:

MM = ComplexF64[1.0237067988305515 -0.03490181882982961 + 0.016873104918059822im 0.018801862132469378 + 0.03132147479849042im; -0.03490181882982961 + 0.016873104918059822im 1.156680064998401 -0.12907274302681904; 0.018801862132469378 + 0.03132147479849042im -0.12907274302681904 0.9866152863493339];
#MM = ComplexF64[1 1im; 1im 5]

show(stdout, "text/plain", MM)
println()
println("Symmetric? ",norm(MM - transpose(MM)) == 0.0)
println("Eigenvalues: ",eigen(MM).values)
LL, _ = my_chol!(copy(MM), LowerTriangular)
show(stdout, "text/plain", LL)
println()
println("Cholesky Error: ",norm(LL * transpose(LL) - MM))

It would be nice to support Symmetric{Complex} matrices.

[ViralBShah]

Doesn't this just need the parent to be extracted and send to the usual cholesky?

[Timeroot]

No, to do a cholesky decomposition on a complex symmetric matrix is a distinct calculation from a cholesky decomposition on a Hermitian matrix. Notably, you can just use +, -, *, / as you would to do a (real symmetric) Cholesky decomposition, but you don't do any complex conjugation. The my_chol! method above is missing the apostrophes (adjoints) from the current version of chol!, and the sqrt(abs(x)) was replaced by just sqrt(x).

[ViralBShah]

Does LAPACK have this (sorry I haven't looked it up)?

[Timeroot]

No, it doesn't. Maybe I should bother them first... :) but it would still be nice to have in Julia too.

[ViralBShah]

Maybe a good place to add that is the GenericLinearAlgebra.jl package to start with. Don't we need some BLAS like stuff to make this fast?

[araujoms]

I don't see what is the problem.

julia> A = Symmetric(randn(ComplexF64, 5, 5));

julia> F = factorize(A); #same thing as bunchkaufman(A)

julia> b = randn(ComplexF64, 5);

julia> x = F \ b;

julia> norm(A * x - b)
1.896870743778145e-15

Julia is doing the same thing as LAPACK here: ZSYSV also solves linear systems of complex symmetric matrices by computing the Bunch-Kaufman decomposition.

Note also that a complex symmetric matrix cannot be positive definite; if a complex matrix is positive definite it is necessarily Hermitian.