Skip to content

Complex element type with SymWoodbury #54

@KlausC

Description

@KlausC

Assuming A and C are symmetric (real or complex) matrices.
In the real case, SymWoodbury(A, U, C) represents A + U * C * U' , because U' == transpose(U).

  1. If U is complex, it is unclear, if A + U * C * U' or A + U * C * transpose(U) is meant; the former matrix is not symmetric.

  2. Wouldn't it be better to stick to the latter definition for the symmetric case and introduce HermWoodbury(A, U, C), and check if Aand C are hermitian?

  3. The checks for symmetry of A and C are not relevant for the validity of the formula, so why are they performed?

  4. If a complex A is not symmetric but hermitian, the checks in SymWoodbury can be easily fooled by providing bunchkaufman(A) or cholesky(A) as the first argument.

  5. I am not sure, if the generalization using pinv is valid in the symmetric but not real case.

From this confusion result the following errors / inconsistencies:

julia> A = Symmetric(rand(2,2) + rand(2,2)*im);

julia> U = rand(2) + rand(2)*im;

julia> C = [1+im;;]

julia> W = SymWoodbury(A, U, C);

julia> Matrix(W) - (A + U * C * U') # that was expected to be zero
2×2 Matrix{ComplexF64}:
       0.0+0.0im       0.246162-0.246162im
 -0.246162+0.246162im       0.0+0.0im

julia> Matrix(W) - (A + U * C * transpose(U)) # that is another candidate, but wasn't either
2×2 Matrix{ComplexF64}:
 0.418788+0.374783im  1.17165+0.582077im
  1.17165+0.582077im  3.13325+0.742344im

julia> inv(W) # bug?
ERROR: MethodError: no method matching -(::Matrix{ComplexF64}, ::ComplexF64)
For element-wise subtraction, use broadcasting with dot syntax: array .- scalar
...

julia> W \ I # that should be the same as inv(W)
2×2 Matrix{ComplexF64}:
 -0.596091-2.90368im  -0.102938+1.71074im
  0.703646+1.45617im  0.0840259-1.23677im

julia> (A + U *C * U') \ I
2×2 Matrix{ComplexF64}:
 -0.596091-2.90368im  -0.102938+1.71074im
  0.703646+1.45617im  0.0840259-1.23677im

julia> (A + U *C * U') \ I - W \ I # surprise! while the matrices differ, their inverses are equal
2×2 Matrix{ComplexF64}:
 -8.88178e-16+1.33227e-15im   5.55112e-16-8.88178e-16im
  7.77156e-16-6.66134e-16im  -3.88578e-16-2.22045e-16im 

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions