I want to find a complex valued matrix A, which fulfills the property A * Transpose(A) = B. The matrix B is complex and symmetric (not hermitian), i.e. B = Transpose(B).
I am working in Julia. I am trying following optimization code which does not work (and is also extremely slow):
using LinearAlgebra
using Optim
# Define the objective function for optimization
function objective(A)
return sum(abs2.(A * transpose(A) - B)) # Sum of squared element wise differences
end
# Initialize a random initial guess for A
n = size(B, 1)
initial_guess = rand(ComplexF64, n, n)
# Perform optimization
result = optimize(objective, initial_guess, BFGS(), Optim.Options(show_trace=true))
# Get the optimized matrix_B_eta from the result
matrix_A_optimized = reshape(result.minimizer, n, n)
How could I constract matrix A?
Here's an edit to include a different, better (i.e. it should be faster and more robust) method based upon the Bunch-Kaufman matrix factorisation. I've left the original post below so the comments still make some kind of sense.
As StackOverflow for some stupid reason doesn't support MathJax it's a bit of a pain to write down the maths of what is going on, but the essence is
Factor B via Bunch-Kaufman: B=Transpose(P)*U*D*Transpose(U)*P
where P is a permutation matrix
U is an upper triangular matrix
D is a block diagonal matrix. The blocks are all 1x1 or 2x2
Then A=Transpose(P)*U*Sqrt(D) we have B=A*Transpose(A)
Note Transpose in the above really is transpose, no complex conjugation is going on.
Here's the code. Again I'll say I am very much a Julia novice, so any ways to improve it would be much appreciated - the Sqrt(D ) function in particular isn't very elegant, but it works and won't be time critical I'm not too worried. That said compiling in the matrix sqrt function does seem to take an inordinate amount of time, especially as it is only used for 2x2 matrices:
using LinearAlgebra
function SymmetricFactorComplexSymmetricMatrix( B::Symmetric{ T } ) where { T <: Complex }
# Main driver
# Given a complex symmetric matrix B return A such that A * A(T)=B
# Method used is via Bunch-Kaufman factorisation
# Factor B such that B=P'UDU'P where
# P is a permutation matrix
# U is a upper triangular matrix
# D is a block diagonal matrix, the blocks being either 1x1 or 2x2
# Due to B being symmetric, D is also symmetric
D, U, p = bunchkaufman( B )
# Form A = P'U*Sqrt(D)
D = SymTridiagonal( D )
A = PermuteResult( U * CalcSqrt2x2BlockDiagonalD( D ), p )
return A
end
function CalcSqrt2x2BlockDiagonalD( D::SymTridiagonal{ T } ) where { T <: Complex }
# Find the square root of the D matrix reurned by the Bunch-Kaufman method
dv = Vector{ T }( undef, size( D, 1 ) )
ev = Vector{ T }( undef, size( D, 1 ) - 1 )
i = 1
while i < size( D, 1 )
# Work out if this is a 2x2 or 1x1 block
D2x2Block = @view D[ i:i + 1, i: i + 1 ]
if abs( D2x2Block[ 1, 2 ] ) ≉ 0.0
# A 2x2 block. Take the sqrt and poke it into the
# arrays that will form the result
SqrtBlock = sqrt( D2x2Block )
dv[ i ] = SqrtBlock[ 1, 1 ]
ev[ i ] = SqrtBlock[ 2, 1 ]
dv[ i + 1 ] = SqrtBlock[ 2, 2 ]
if i ≠ size( D, 1 ) - 1
ev[ i + 1 ] = 0.0
end
i = i + 2
else
# A 1x1 block. Take the sqrt and poke it into the
# arrays that will form the result
dv[ i ] = sqrt( D[ i, i ] )
ev[ i ] = 0.0
i = i + 1
end
end
# If required catch the last element. Must be a 1x1 block
if i == size( D, 1 )
dv[ i ] = sqrt( D[ i, i ] )
end
return SymTridiagonal( dv, ev )
end
function PermuteResult( Fin::Matrix{ T }, p::Vector{ IT } ) where { T <: Complex, IT <: Integer }
# Permute the Fin matrix accordin to the permutation given by p
Fout = similar( Fin )
for i in axes( Fin, 2 )
Fout[ p, i ] = Fin[ :, i ]
end
return Fout
end
# Checking functions
function CheckIt( A::Symmetric{ T } ) where{ T <: Complex }
F = SymmetricFactorComplexSymmetricMatrix( A )
return maximum( abs.( A - F * ComplexTranspose( F ) ) )
end
function ComplexTranspose( Q::Matrix )
QT = similar( Q )
for j in axes( Q, 1 )
for i in axes( Q, 2 )
QT[ i, j ] = Q[ j, i ]
end
end
return QT
end
function ManyTest( ntests::Integer, maxsize = 1000::Integer )
maxerror = 0.0
for i = 1:ntests
n = Int( floor( maxsize * rand() + 1 ) )
A = Symmetric( rand( ComplexF64, n, n ) )
error = CheckIt( A )
println( "Test ", i, " size ", n, " error ", error )
maxerror = max( error, maxerror )
end
printstyled( "Max error: ", maxerror, "\n"; color = :lightred )
return maxerror
end
And here's it being tested
ijb@LAPTOP-GUG8KQ9I:~/julia/test$ julia -O3
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.9.0 (2023-05-07)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> using Revise, OhMyREPL, Debugger, Test, Pkg
julia> include( "CompSym2.jl")
ManyTest (generic function with 2 methods)
julia> ManyTest(25,5000)
Test 1 size 3411 error 8.1594242764668e-13
Test 2 size 1484 error 5.172734121734416e-13
Test 3 size 4184 error 8.575452117894415e-13
Test 4 size 2018 error 4.830629499161676e-13
Test 5 size 1658 error 6.091354345513839e-13
Test 6 size 2841 error 1.1587510250913577e-12
Test 7 size 1589 error 7.092154735016015e-13
Test 8 size 4276 error 1.2948552601956827e-12
Test 9 size 3435 error 1.0944800551535788e-12
Test 10 size 1075 error 5.827720133625812e-13
Test 11 size 3237 error 1.1542970580509046e-12
Test 12 size 861 error 1.9065513187377532e-13
Test 13 size 946 error 2.5678632958931523e-13
Test 14 size 3953 error 1.2280159009015767e-12
Test 15 size 4737 error 8.789983456404833e-13
Test 16 size 4477 error 2.315395050493212e-12
Test 17 size 3024 error 5.390843920673155e-13
Test 18 size 3610 error 1.9573189511434918e-12
Test 19 size 273 error 9.825630577843e-14
Test 20 size 1170 error 6.206737441326594e-13
Test 21 size 1511 error 5.482164493856334e-13
Test 22 size 770 error 2.495593681148973e-13
Test 23 size 738 error 2.9098725209054565e-13
Test 24 size 2692 error 1.059692174415002e-12
Test 25 size 2288 error 5.528934962642065e-13
Max error: 2.315395050493212e-12
2.315395050493212e-12
Original Method:
Here's an idea which seems to work but note the caveat below. Also note I am very much a Julia novice ( but very experienced in Fortran and C ), so this might be a good example of being able to write Fortran in any language ...
function ComplexTranspose( Q::Matrix )
QT = similar( Q )
for j in axes( Q, 1 )
for i in axes( Q, 2 )
QT[ i, j ] = Q[ j, i ]
end
end
return QT
end
function GenQOrtho( Q::Matrix )
QOrtho = copy( Q )
for i in axes( QOrtho, 2 )
α = sum( QOrtho[ :, i ] .* QOrtho[ :, i ] )
α = sqrt( α )
QOrtho[ :, i ] = QOrtho[ :, i ] .* ( 1.0 / α )
end
return QOrtho
end
function SymmetricFactorComplexSymmetricMatrix( A::Matrix{ T } ) where { T <: Complex }
λ, Q = eigen( A )
QOrtho = GenQOrtho( Q )
return Diagonal( sqrt.( λ ) ) * ComplexTranspose( QOrtho )
end
function CheckIt( A::Matrix )
F = SymmetricFactorComplexSymmetricMatrix( A )
return maximum( abs.( A - ComplexTranspose( F ) * F ) )
end
function ManyTest( ntests::Integer )
maxerror = 0.0
for i = 1:ntests
n = Int( floor( 100 * rand() + 1 ) )
A = Matrix( Symmetric( rand( ComplexF64, n, n ) ) )
error = CheckIt( A )
println( "Test ", i, " size ", n, " error ", error )
maxerror = max( error, maxerror )
end
printstyled( "Max error: ", maxerror, "\n"; color = :lightred )
return maxerror
end
Running in the REPL:
ijb@LAPTOP-GUG8KQ9I:~/julia/test$ julia
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.9.0 (2023-05-07)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> using Revise, OhMyREPL, Debugger, Test, Pkg, LinearAlgebra
julia> include( "CompSym.jl")
ManyTest (generic function with 1 method)
julia> ManyTest(25)
Test 1 size 83 error 9.944536531374647e-14
Test 2 size 11 error 4.491946442396235e-15
Test 3 size 4 error 8.95090418262362e-16
Test 4 size 48 error 5.541685356342546e-14
Test 5 size 58 error 6.776845358219464e-14
Test 6 size 83 error 1.2052768509123334e-13
Test 7 size 23 error 1.602186845968893e-14
Test 8 size 79 error 1.1038213408936726e-13
Test 9 size 77 error 6.335118026431236e-13
Test 10 size 27 error 1.5122688478052736e-14
Test 11 size 62 error 1.3015134857799234e-13
Test 12 size 59 error 9.97847338085134e-14
Test 13 size 36 error 3.466048004558062e-14
Test 14 size 28 error 1.9942028578972164e-14
Test 15 size 74 error 1.899255854395612e-13
Test 16 size 13 error 4.228354934573259e-15
Test 17 size 83 error 3.63646988782893e-13
Test 18 size 56 error 6.795807448681065e-14
Test 19 size 89 error 1.6855954444400822e-13
Test 20 size 97 error 8.38314879211208e-14
Test 21 size 69 error 8.592074666261367e-14
Test 22 size 7 error 3.5283439685664285e-15
Test 23 size 5 error 1.6653345369377348e-15
Test 24 size 60 error 6.41406310718282e-14
Test 25 size 35 error 2.6668010011248128e-14
Max error: 6.335118026431236e-13
It's based on the fact that the eigenvectors of a complex symmetric matrix are orthogonal in a transpose sense and not in the usual conjugate transpose sense. Thus we can use the normal eigenvalue machinery, but have to carefully renormalise the eigenvectors so that the length of the vectors are 1 in the transpose as opposed to the adjoint sense - which is a bit of a pain as most of the Julia tools such as dot() or norm() assume you want the conjugate (as you usually do), so here you have to do some work "manually". I suspect this is the reason your solution doesn't work - the transpose function is performing a conjugation as it sees a complex matrix.
The caveat is I am not 100% convinced this will always work when you have degeneracies. The problem is Julia doesn't seem to have a specialised Complex Symmetric diagonaliser, and so I can't see why evecs corresponding to degenerate evals are guaranteed to be orthogonal (in any sense). I can see how to fix this (diagonalise in the subspace), but currently I can't generate an example of this, and am too tired to write the code anyway. I'll see if I can over the weekend.