# What is the point of the Symmetric type?

**URL:** <https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773>\
**Category:** Internals & Design\
**Tags:** linearalgebra\
**Created:** [June 4, 2020, 11:21pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773 "2020-06-04T23:21:17Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Grayscale](https://avatars.discourse-cdn.com/v4/letter/g/e19b73/32.png) [@Grayscale](https://discourse.julialang.org/u/Grayscale)\
**Post date:** [June 4, 2020, 11:21pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/1 "2020-06-04T23:21:18Z")

</div>

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](https://stackoverflow.com/questions/62058109/what-is-the-point-of-the-symmetric-type-in-julia)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 5, 2020, 1:31am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/2 "2020-06-05T01:31:28Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Grayscale](https://avatars.discourse-cdn.com/v4/letter/g/e19b73/32.png) [@Grayscale](https://discourse.julialang.org/u/Grayscale)\
**Post date:** [June 5, 2020, 1:35am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/3 "2020-06-05T01:35:20Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [June 5, 2020, 1:42am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/4 "2020-06-05T01:42:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 5, 2020, 1:47am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/5 "2020-06-05T01:47:48Z")

</div>

> [@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](http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=3592#p9189). 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.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [June 5, 2020, 3:06am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/6 "2020-06-05T03:06:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [June 5, 2020, 8:20am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/7 "2020-06-05T08:20:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 5, 2020, 10:16am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/8 "2020-06-05T10:16:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 5, 2020, 12:03pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/9 "2020-06-05T12:03:51Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![CodeLenz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/codelenz/32/3419_2.png) [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)\
**Post date:** [June 5, 2020, 1:21pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/10 "2020-06-05T13:21:34Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![platawiec](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/platawiec/32/31914_2.png) [@platawiec](https://discourse.julialang.org/u/platawiec)\
**Post date:** [June 5, 2020, 2:26pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/11 "2020-06-05T14:26:17Z")

</div>

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](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.

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 5, 2020, 2:44pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/12 "2020-06-05T14:44:25Z")

</div>

> [@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_](https://en.wikipedia.org/wiki/Occam%27s_razor).

We work around this in `LinearAlgebra` by having aliases like:

```julia
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.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 5, 2020, 3:55pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/13 "2020-06-05T15:55:34Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![CodeLenz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/codelenz/32/3419_2.png) [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)\
**Post date:** [June 5, 2020, 4:38pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/14 "2020-06-05T16:38:45Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [June 5, 2020, 9:43pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/16 "2020-06-05T21:43:26Z")

</div>

> [@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:

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

```

together with

```julia
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.

---

<div class="post-metadata">

**Author:** ![Grayscale](https://avatars.discourse-cdn.com/v4/letter/g/e19b73/32.png) [@Grayscale](https://discourse.julialang.org/u/Grayscale)\
**Post date:** [June 10, 2020, 5:17am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/17 "2020-06-10T05:17:13Z")

</div>

By the way, in looking into the [LAPACK documentation](http://www.icl.utk.edu/~mgates3/docs/lapack.html), 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).

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [June 13, 2020, 3:41am UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/18 "2020-06-13T03:41:28Z")

</div>

> [@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):

```julia
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:

```julia
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

```

```julia
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:

```julia
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:  
 ![solution_time](https://global.discourse-cdn.com/julialang/original/3X/e/2/e2ebe50d29ed34c6891b1865d2cb22cddc3f6c4f.png)  
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.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [April 26, 2023, 8:20pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/19 "2023-04-26T20:20:07Z")

</div>

> [@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](http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=3592#p9189). 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?

---

<div class="post-metadata">

**Author:** ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)\
**Post date:** [April 26, 2023, 8:34pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/20 "2023-04-26T20:34:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [April 26, 2023, 9:37pm UTC](https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/21 "2023-04-26T21:37:56Z")

</div>

Thanks!