# What is the point of the Symmetric type?

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

1 Like

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.

11 Likes

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…

1 Like

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.

1 Like

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.

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.

5 Likes

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.

3 Likes

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.

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.
1 Like

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.

2 Likes

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.

3 Likes

There are some variants of look-ahead Lanczos decomposition algorithms which have simplifications in the complex-symmetric case , 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: https://web.cs.ucdavis.edu/~bai/ET/lanczos_methods/overview_QMRPACK.html

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

 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. https://doi.org/10.1137/0914009
 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. https://doi.org/10.1137/0915022
 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. https://doi.org/10.1137/0913023

2 Likes

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.

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.

6 Likes

@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.

1 Like

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

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.

4 Likes

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).

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
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))
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.

1 Like