What is the point of the Symmetric type?

[Grayscale]

Cross-posting this question at the suggestion of a commenter:

https://stackoverflow.com/questions/62058109/what-is-the-point-of-the-symmetric-type-in-julia

[stevengj]

Always use Hermitian. For real elements, there is no penalty compared to Symmetric.

There aren’t any specialized routines for complex Symmetric matrices that I know of. My feeling is that it was probably a mistake to have a separate Symmetric type in LinearAlgebra, but it is hard to remove at this point.

[Grayscale]

Hm interesting. Thanks!

Do you know why then there is a SymTridiagonal type but no HermTridiagonal? It seems like the latter would be more useful…

[Oscar_Smith]

I’d argue that Symmetric is important for the significant portion of people who use LinearAlgebra but don’t do anything with complex numbers or have a significant Mathy background.

[stevengj]

@Grayscale

Do you know why then there is a SymTridiagonal type but no HermTridiagonal ?

I tend to agree that we should have HermTridiagonal instead, but you should realize that any Hermitian tridiagonal matrix T can be converted to a similar real-symmetric tridiagonal matrix T' = D^* T D by a unitary diagonal scaling matrix D — see e.g. the discussion here. So, in principle, you can always work with real tridiagonal matrices.

For the same reason LAPACK’s eigensolver for complex Hermitian matrices first reduce them to similar real-symmetric tridiagonal matrices (this is done by the hessenberg factorization routine in Julia), and it was to reflect these routines that SymTridiagonal was first introduced.

@Oscar_Smith

I’d argue that Symmetric is important for the significant portion of people who use LinearAlgebra but don’t do anything with complex numbers or have a significant Mathy background.

Those people can just use Hermitian too. For real matrices, Hermitian and Symmetric are equivalent. Even if they don’t know the word “Hermitian,” learning a new word is not such a big obstacle.

[PeterSimon]

From your SO post:

it would be surprising to want to have a symmetric as opposed to Hermitian complex matrix (maybe I am mistaken on that though).

Complex symmetric matrices occur in electrical engineering applications. In circuit theory, multiport Impedance matrices, admittance matrices, and scattering matrices of reciprocal networks are all complex symmetric. LAPACK routine zsysv solves a system of linear equations with a complex symmetric coefficient matrix.

[simeonschaub]

Are there functions that actually see a significant speedup with the Symmetric/Hermitian wrapper? The functions I tried didn’t really get any faster and matrix-vector multiply got a lot slower, for example.

[Tamas_Papp]

@PeterSimon

Complex symmetric matrices occur in electrical engineering applications.

All kinds of specialized matrices exist for various applications, the important question is whether they should live in LinearAlgebra. In this case, arguably they should not, since 99% of users don’t need the distinction beyond Hermitian and may just find it confusing.

It would be better to

1. deprecate LinearAlgebra.Symmetric in the long run (when we get a chance), then remove it altogether,
2. have a package for complex symmetric matrices and their specialized methods.

[stevengj]

@PeterSimon

Complex symmetric matrices occur in electrical engineering applications.

Yes, Lorentz reciprocity makes matrices that arise from Maxwell’s equations complex-symmetric. However, since we have no algorithms that take advantage of the complex-symmetric case, there is no benefit to using Symmetric here.

[CodeLenz]

I use it to enforce symmetry when performing Cholesky factorization. My matrices are (theoretically) symmetric but very small round errors are detected as non symmetry by cholesky(A). So cholesky(Symmetric(A)) does the job.

[platawiec]

There are some variants of look-ahead Lanczos decomposition algorithms which have simplifications in the complex-symmetric case [1], at least that’s how Freund and Nachtigal implemented them. QMR and other iterative solver methods were developed for solving this class of systems [2, 3]. See qmrpack: QMRPACK

Some textbook implementations don’t follow their algorithm, however. It is a somewhat niche method as well.

[1] Freund, R. W., Gutknecht, M. H., & Nachtigal, N. M. (1993). An Implementation of the Look-Ahead Lanczos Algorithm for Non-Hermitian Matrices. SIAM Journal on Scientific Computing , 14 (1), 137–158. The Society for Industrial and Applied Mathematics
[2] Freund, R. W., & Nachtigal, N. M. (1994). An Implementation of the QMR Method Based on Coupled Two-Term Recurrences. SIAM Journal on Scientific Computing , 15 (2), 313–337. The Society for Industrial and Applied Mathematics
[3] Freund, R. W. (1992). Conjugate Gradient-Type Methods for Linear Systems with Complex Symmetric Coefficient Matrices. SIAM Journal on Scientific and Statistical Computing , 13 (1), 425–448. The Society for Industrial and Applied Mathematics

[stevengj]

@CodeLenz

I use it to enforce symmetry when performing Cholesky factorization.

Cholesky factorization only works for Hermitian matrices. So you can do cholesky(Hermitian(A)) in all cases; this is equivalent to cholesky(Symmetric(A)) for real A, while for complex A doing cholesky(Symmetric(A)) is incorrect.

@platawiec

There are some variants of look-ahead Lanczos decomposition algorithms

Yes, I know there are specialized algorithms for complex-symmetric matrices in some cases. But they don’t exist in LAPACK or in the LinearAlgebra package. If someone were to implement a package of specialized algorithms for complex-symmetric matrices, that would be a good place for a Symmetric type to live.

The problem with having it in LinearAlgebra, with no algorithms that exploit it, is that it just causes confusion to have two types (Hermitian and Symmetric) that are equivalent for the real case and only one of them is usable in the complex case. Entia non multiplicanda sunt.

We work around this in LinearAlgebra by having aliases like:

const RealHermSymComplexHerm{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}

but it would be nicer to just be able to say Hermitian.

At this point it is water under the bridge, however — it is pretty clear that we can’t move Symmetric into a package until Julia 2.0, and it might be too disruptive to do it even then. It is probably more productive to talk about implementing a package that would actually make Symmetric{<:Complex} matrices useful.

My main point is to answer the original poster: there is zero difference in practice between Symmetric and Hermitian for real matrices, and in fact they share the same code under the hood. Probably this should be clarified in the documentation.

[Elrod]

@CodeLenz’s example is also a case where the wrapper provides a substantial performance benefit, via skipping the check (although using Hermitian instead of Symmetric, per the above suggestion):

julia> using LinearAlgebra, BenchmarkTools

julia> S = rand(2000,3000) |> x -> x * x';

julia> hS = Hermitian(S);

julia> @benchmark cholesky($S)
BenchmarkTools.Trial:
  memory estimate:  30.52 MiB
  allocs estimate:  6
  --------------
  minimum time:     13.474 ms (0.00% GC)
  median time:      13.995 ms (0.00% GC)
  mean time:        15.298 ms (3.95% GC)
  maximum time:     79.703 ms (76.94% GC)
  --------------
  samples:          327
  evals/sample:     1

julia> @benchmark cholesky($hS)
BenchmarkTools.Trial:
  memory estimate:  30.52 MiB
  allocs estimate:  6
  --------------
  minimum time:     7.709 ms (0.00% GC)
  median time:      8.282 ms (0.00% GC)
  mean time:        9.273 ms (4.87% GC)
  maximum time:     69.711 ms (88.67% GC)
  --------------
  samples:          539
  evals/sample:     1

I think they’re also convenient for dispatch.
For example, you may define

multivariate_normal_logpdf(y, μ, Σ::Hermitian) = # covariance matrix argument
multivariate_normal_logpdf(y, μ, U::UpperTriangular) = # cholesky factor of covariance matrix argument
multivariate_normal_logpdf(y, μ, L::LowerTriangular) = # cholesky factor of covariance matrix argument

You may also see gains from not needing to calculate the entire matrix.
Perhaps the entries of the matrix are defined by a function such as f(x,i,j) = exp(-abs(x[i]-x[j])); it’s convenient not to have to manually copy one triangle to the other.

[CodeLenz]

I am not advocating for it, just providing an user case for the Symmetric/Hermitian type.

[juthohaegeman]

@stevengj

it is pretty clear that we can’t move Symmetric into a package until Julia 2.0, and it might be too disruptive to do it even then. It is probably more productive to talk about implementing a package that would actually make Symmetric{<:Complex} matrices useful.

Could one not do:

const Symmetric{T<:Real, A<:AbstractMatrix{T}} = Hermitian{T,A}

together with

Symmetric(a::AbstractMatrix{<:Real}) = Hermitian(a)
Symmetric(a::AbstractMatrix{<:Complex}) = ... (some deprecation warning and/or referral to a package that provides actual support for ComplexSymmetric matrices)

It seems like this would not be breaking at all for most use cases, which simply involve wrapping a matrix. That type alias could stay one indefinitely rather than only as a transition period, for those users who find Hermitian to math-like.

[Grayscale]

By the way, in looking into the LAPACK documentation, I noticed both that complex symmetric matrices come up somewhat often and that symmetric tridiagonal matrices, but not Hermitian tridiagonal matrices, are featured. This is me speculating, but it thus seems plausible that the Julia LinearAlgebra package has in mind the LAPACK types in its design (even though it appears not to currently use any of the associated routines for e.g. complex symmetric matrices).

[PeterSimon]

@stevengj

… there are specialized algorithms for complex-symmetric matrices in some cases. But they don’t exist in LAPACK …

I was curious after reading this if @stevengl knew some reason for not using zsysv in LAPACK. I created a wrapper in Julia copying the examples in LinearAlgebra.LAPACK (the following Julia code blocks are copied from a Jupyter notebook):

import LinearAlgebra.LAPACK: chkuplo, chkargsok, chknonsingular, BlasInt, chkstride1, checksquare

function zsysv!(uplo::AbstractChar, A::AbstractMatrix{ComplexF64}, B::AbstractVecOrMat{ComplexF64})
    Base.require_one_based_indexing(A, B)
    chkstride1(A,B)
    n = checksquare(A)
    chkuplo(uplo)
    if n != size(B,1)
        throw(DimensionMismatch("B has first dimension $(size(B,1)), but needs $n"))
    end
    ipiv  = similar(A, BlasInt, n)
    work  = Vector{ComplexF64}(undef, 1)
    lwork = BlasInt(-1)
    info  = Ref{BlasInt}()
    for i = 1:2  # first call returns lwork as work[1]
        ccall((:zsysv64_, Base.liblapack_name), Cvoid,
            (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, Ptr{BlasInt},
            Ptr{ComplexF64}, Ref{BlasInt}, Ptr{ComplexF64}, Ref{BlasInt}, Ptr{BlasInt}),
            uplo, n, size(B,2), A, max(1,stride(A,2)), ipiv, B, max(1,stride(B,2)),
            work, lwork, info)
        chkargsok(info[])
        chknonsingular(info[])
        if i == 1
            lwork = BlasInt(real(work[1]))
            resize!(work, lwork)
        end
    end
    B, A, ipiv
    return
end

A quick test:

import LinearAlgebra: norm
n = 2000
A0 = complex.(rand(n,n), rand(n,n))
A0 = A0 + transpose(A0) # Force A to be complex-symmetric
B0 = complex.(rand(n), rand(n))
A = copy(A0); B = copy(B0)
zsysv!('U', A, B)
norm(A0*B - B0)  #  Result is 7.283816110456216e-11

x = A0\B0;
norm(A0*x - B0)  # Result is 3.566949817832705e-11

so the implementation seems to be roughly as accurate as the algorithm selected by \. Now for some timing runs:

using BenchmarkTools
orders = 10*(2 .^(0:8))
t_backslash = zeros(length(orders))
t_zsysv = zeros(length(orders))
for (i,n) ∈ enumerate(orders)
    A0 = complex.(rand(n,n), rand(n,n))
    A0 = A0 + transpose(A0) # Force A to be complex-symmetric
    B0 = complex.(rand(n), rand(n))
    t_backslash[i] = @belapsed $A0 \ $B0
    t_zsysv[i] = @belapsed zsysv!('U', A, B) setup=(A = copy($A0); B = copy($B0))
end

and after plotting the results:

I was very surprised to see that for matrix orders above 320, zsysv! is 2 to 3 times slower than the general complex matrix algorithm selected by the backslash operator. Did I screw up the benchmarking or the implementation for zsysv!? Because I expected that the specialized algorithm in LAPACK for symmetric matrices would be faster than that for general complex matrices, for all matrix orders.

[Mason]

@stevengj

I tend to agree that we should have HermTridiagonal instead, but you should realize that any Hermitian tridiagonal matrix TTT can be converted to a similar real-symmetric tridiagonal matrix T’ = D^* T DT′=D∗TDT’ = D^* T D by a unitary diagonal scaling matrix DDD — see e.g. the discussion here . So, in principle, you can always work with real tridiagonal matrices.

That link is dead now, are you aware of other useful sources on this?

[jar1]

That link is archived
https://web.archive.org/web/20201028202318/http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=3592#p9189

[Mason]

Thanks!