Complex diagonal elements of APA'

[nan2ge1]

Suppose P is a real square diagonal matrix, then APA' should be Hermitian with real diagonal elements. There is however a chance for some diagonal elements to have nonzero imaginary parts.

julia> using LinearAlgebra, SparseArrays

julia> m, n = 2, 3;

julia> A = randn(m,n) + im*randn(m,n);

julia> p = abs.( randn(n) );

julia> A * spdiagm(p) * A'
2×2 Array{Complex{Float64},2}:
 1.34723+5.55112e-17im               -0.18145+1.30975im
         -0.18145-1.30975im  4.48512+2.22045e-16im     

Enforcing Hermitian (e.g. as input for cholfact) will end up with the following error.

julia> Hermitian( A * spdiagm(p) * A' )
ERROR: ArgumentError: Cannot construct Hermitian from matrix with nonreal diagonals
 in Hermitian{T,S<:AbstractArray{T,2}}(::Array{Complex{Float64},2}, ::Symbol) at ./linalg/symmetric.jl:50
 in Hermitian{T,S<:AbstractArray{T,2}}(::Array{Complex{Float64},2}) at ./linalg/symmetric.jl:48

A work-around could be to re-assign the diagonal elements to their real parts, but that would also slow things down.

julia> versioninfo()
Julia Version 0.5.0
Commit 3c9d753 (2016-09-19 18:14 UTC)
Platform Info:
  System: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
  WORD_SIZE: 64
  BLAS: libopenblas (NO_LAPACKE DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: liblapack.so.3
  LIBM: libopenlibm
  LLVM: libLLVM-3.7.1 (ORCJIT, broadwell)

[StefanKarpinski]

Is the general "contract" of these matrix wrappers to leave the underlying data untouched?

[andreasnoack]

This is indeed a problem. See JuliaLang/julia#17367

[nan2ge1]

@StefanKarpinski FYI: it seems that MATLAB does no "post-processing".

>> m = 2;
>> n = 3;
>> A = randn(m,n) + 1j*randn(m,n);
>> p = abs( randn(n,1) );
>> A * diag(p) * A'
ans =
   2.8503 - 0.0000i  -2.1165 + 0.4488i
  -2.1165 - 0.4488i  10.8608 - 0.0000i
>> ishermitian( A * diag(p) * A' )
ans =
  logical
   0
>> chol( A * diag(p) * A' )
Error using chol
Matrix must be positive definite with real diagonal. 

[nan2ge1]

@andreasnoack Thanks very much for the cross-reference. I like the idea of using (A+A')/2 for Hermitian A.

[dkarrasch]

This was closed by JuliaLang/julia#17367.