Add matrix symmetrizing functions

[ararslan]

This adds functions symmetrize and symmetrize!, which produce symmetric matrices based on the input X via efficient computation of (X + X') / 2.

Benchmark comparison with a naive implementation:

julia> @benchmark symmetrize!(X) setup=(X = randn(20, 20))
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     120.718 ns (0.00% GC)
  median time:      121.681 ns (0.00% GC)
  mean time:        124.094 ns (0.00% GC)
  maximum time:     263.909 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     905

julia> f!(X) = (X .= (X .+ X') ./ 2; X)
f! (generic function with 1 method)

julia> @benchmark f!(X) setup=(X = randn(20, 20))
BenchmarkTools.Trial: 
  memory estimate:  3.25 KiB
  allocs estimate:  1
  --------------
  minimum time:     539.921 ns (0.00% GC)
  median time:      587.149 ns (0.00% GC)
  mean time:        668.309 ns (10.85% GC)
  maximum time:     152.017 μs (99.49% GC)
  --------------
  samples:          10000
  evals/sample:     191

Currently this only deals with real matrices. I can add complex as well if that'd be useful (assuming this is something we want in general).

[jebej]

A complex hermitianize! would definitely be useful!

[andreasnoack]

A complex hermitianize! would definitely be useful!

I think we should stretch the terminology a bit here and have a flag in the symmetrize function for conjugation which should default to true. It will, by far, be the more common operation.

[fredrikekre]

Can we call this symmetric instead?

[ararslan]

That was my first instinct but we already have such a function, it just isn't exported.

help?> LinearAlgebra.symmetric
  symmetric(A, uplo=:U)

  Construct a symmetric view of A. If A is a matrix, uplo controls whether the upper (if uplo = :U) or lower (if uplo = :L) triangle of A is used to
  implicitly fill the other one. If A is a Number, it is returned as is.

  If a symmetric view of a matrix is to be constructed of which the elements are neither matrices nor numbers, an appropriate method of symmetric has to be
  implemented. In that case, symmetric_type has to be implemented, too.

[fredrikekre]

Okay, can we rename that internal function then? We use symmetric in https://github.com/KristofferC/Tensors.jl for this purpose for example.

[ararslan]

Okay, can we rename that internal function then?

We'd need to do a PkgEval run to determine whether it's being used, as that has the potential to be breaking.

[ararslan]

Looking at it more, it would be kind of weird to change the internal one, as it's named in keeping with the convention of lowercase/uppercase symmetry. That is, symmetric constructs a Symmetric (and likewise for hermitian) but also accepts numbers and is overloadable for other types.

[musm]

shouldn't these function then be returning Symmetric types after 'symmetrizing' them

[ararslan]

IMO no, because it doesn't make sense for symmetrize!(X) to not return X, and it would be weird for symmetrize to not have the same behavior as its in-place counterpart.

[andreasnoack]

it doesn't make sense for symmetrize!(X) to not return X

I don't agree. E.g. cholesky!(X) doesn't return X. However, I'm not completely sure what I prefer. I think I'm slightly in favor of returning a Symmmetric/Hermitian. Then we could also save (slightly more than) half the stores.

[ararslan]

I made the change locally but it looks like the Symmetric and Hermitian constructors require 1-based indexing:

  Got exception outside of a @test
  AssertionError: !(has_offset_axes(data))
  Stacktrace:
   [1] Type at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/symmetric.jl:9 [inlined]
   [2] Type at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/symmetric.jl:50 [inlined]
   [3] Type at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/symmetric.jl:49 [inlined]
   [4] #symmetrize!#3 at ./REPL[2]:13 [inlined]

As it stands, symmetrize! explicitly allows offset axes.

[stevengj]

Can we just call it hermitian(A) and hermitian!(A) and remove the non-conjugated option?

(In hindsight, I wish we'd just gotten rid of Symmetric completely — for the real case it's equivalent to Hermitian and forces us to HermOrSym everywhere, and for the complex case it's nearly useless. Maybe we can at least stop adding new [Ss]ymmetric functions.)

[ararslan]

I'd be fine with that. Would hermitian always return a Hermitian-wrapped matrix? Or do we want to support offset indices?

[stevengj]

👍 on always returning a Hermitian wrapper. Supporting offset arrays could come later, if ever — most of the linear-algebra routines don't support it, so I'm not sure how worthwhile it is to support offset arrays piecemeal.

[ararslan]

Can we just call it hermitian(A) and hermitian!(A) and remove the non-conjugated option?

Heh, turns out that one exists too. Like symmetric, it's just not exported.

help?> LinearAlgebra.hermitian
  hermitian(A, uplo=:U)

  Construct a hermitian view of A. If A is a matrix, uplo controls whether the upper (if uplo = :U) or lower (if uplo = :L) triangle of A is used to
  implicitly fill the other one. If A is a Number, its real part is returned converted back to the input type.

  If a hermitian view of a matrix is to be constructed of which the elements are neither matrices nor numbers, an appropriate method of hermitian has to be
  implemented. In that case, hermitian_type has to be implemented, too.

[ararslan]

Bump. We can't call it symmetric or hermitian since those functions already exist, so is the current name acceptable? If so, where did we land with returning a wrapper? Yes and ignore non-1-based indices?

[jebej]

I'd be fine with not returning a wrapper, espcially given the in-place function argument above.

Not sure what to do the about the hermitian and symmetric functions, they seem pretty useless. If we want to leave them as is I'm okay with hermitianize or symmetrize. Note that choosing hermitianize might cause discoverability issues with people who work only with real matrices might not know the term "Hermitian".

[ararslan]

Bump. Any interest in this or should it be closed?

[stevengj]

Since LinearAlgebra.hermitian is not exported, aren't we free to rename it to hermitianview or something? Is it documented in the manual?

I'm not a huge fan of the name symmetrize for something that take the Hermitian part, but I definitely think we should have some kind of functions with this functionality in LinearAlgebra.

We could call it hermitianpart, I guess, but I like the name hermitian(a) in analogy to real(a), especially if it returns a Hermitian wrapper which I think would be advisable.

[ararslan]

Is it documented in the manual?

Doesn't appear to be.

Here's a good approximation of all of the uses in LinearAlgebra plus registered packages: https://juliahub.com/ui/Code?q=%5Cbhermitian%5C%28&w=true&r=true. Only MutableArithmetics seems to be using the one from LinearAlgebra; LinearOperators is defining its own to mean something different, which it exports.

[oscardssmith]

This looks useful. Can it be rebased? Added triage label to hopefully give it some attention and so we can decide on names and whether to add.

[stevengj]

Is it possible to have s != 1? a[begin:end] only works if the step is 1, so I thought we required consecutive ranges of indices.

[dkarrasch]

• Was it timed how much is saved by only updating the upper (or lower) part of the matrix, because it will be wrapped in Hermitian? I could see that people want this functionality where the whole matrix is updated, for applications different than LAPACK calls and the Hermitian wrappers. If the gain is negligible, it would be a missed opportunity not to make the full matrix symmetric/hermitian.

This could be keyword argument: full that is false by default. But if it's true, it could be like full && copytri!(...). To make it up to the user?

[Jutho]

This could be keyword argument: full that is false by default. But if it's true, it could be like full && copytri!(...). To make it up to the user?

The keyword is fine, but that implementation seems suboptimal. You would be running over the data (which is already not very cache friendly due to the fact that you need A[i,j] and A[j,i] at the same time) a second time.

I checked the overhead of immediately putting the result in A[j,i] as well. For smaller matrices (100 x 100) the overhead is actually quite significant (contrary to my expectations), but for larger matrices, it goes down.

[stevengj]

Suggested change

	_hermitianpart(a::T) where {T<:Number} = convert(Base.promote_op(/, T, Int), real(a))
	_hermitianpart(a::Number) = real(a)

As far as I can tell there is no reason to do a conversion here.

[ararslan]

I suggested doing so to ensure that typeof(hermitianpart(x)) == eltype(hermitianpart([x;;])) for x::Number

[dkarrasch]

When this is used along the diagonal of the matrix, then that conversion happens anyway, since the target array has that specific eltype. In that case, the diagonal would be filled with complex numbers with vanishing real part, not with real numbers! Ananlogously, computing (a+a') == (a+conj(a))/2 for complex a yields a complex number, in contrast to real(a). OTOH, there is the "operator real part" interpretation, which would perfectly allow the pure real(a) implementation. I believe, however, that we should remove the hermitianpart(a::Number) method, because it is (a) currently undocumented, and (b) never called internally (we call _hermitianpart instead).

[dkarrasch]

I removed hermitianpart(::Number) and the type conversion altogether, because the goal is to provide a "matrix symmetrizing function". hermitianpart(::Number) does not interfere with that goal, so can be an independent PR if necessary. People may then also have a concrete use case and expectations regarding the return type. It would also require an explicit comment in the docstring on how hermitianpart is working on Numbers.

[dkarrasch]

Ready to go?

[thchr]

Would there be any appetite for a PR adding the analogue antihermitianpart(!) along the lines of the opimag described above? I had a need for it now and found the asymmetry unfortunate (mainly for the lack of an in-place version in practice).