我想找到一个复值矩阵 A，它满足属性 A * Transpose(A) = B。矩阵 B 是复数且对称的(不是 Hermitian 矩阵)，即 B = Transpose(B)。
我在 Julia 工作。我正在尝试以下优化代码，但它不起作用(而且速度也非常慢):
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)
如何构造矩阵A？
最佳答案
这里的编辑包括一个不同的、更好的(即它应该更快、更稳健)基于 Bunch-Kaufman 矩阵分解的方法。我已将原始帖子留在下面，因此评论仍然有意义。
由于 StackOverflow 由于某些愚蠢的原因不支持 MathJax，因此写下正在发生的数学运算有点痛苦，但本质是
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.
这是代码。我再说一次，我是一个 Julia 新手，所以任何改进它的方法都将非常感激 - 特别是 Sqrt(D ) 函数不是很优雅，但它可以工作，并且不会对时间至关重要我'我不太担心。也就是说，在矩阵 sqrt 函数中进行编译似乎确实花费了过多的时间，特别是因为它仅用于 2x2 矩阵:
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
这是正在测试的
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
原始方法:
这是一个似乎可行的想法，但请注意下面的警告。另请注意，我是 Julia 新手(但在 Fortran 和 C 方面非常有经验)，所以这可能是能够用任何语言编写 Fortran 的一个很好的例子......
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
结束
在 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
它基于复对称矩阵的特征向量是 orthogonal in a transpose sense and not in the usual conjugate transpose sense 的事实。因此，我们可以使用正常的特征值机制，但必须仔细地重新规范化特征向量，以便向量的长度在转置中为 1，而不是伴随意义 - 这有点痛苦，因为大多数 Julia 工具，例如dot() 或 norm() 假设您想要共轭(就像您通常做的那样)，所以这里您必须“手动”做一些工作。我怀疑这就是您的解决方案不起作用的原因 - transpose 函数正在执行共轭，因为它看到一个复杂的矩阵。
需要注意的是，我并不 100% 相信当你有退化时这总是有效。问题是 Julia 似乎没有专门的复杂对称对角线器，所以我不明白为什么对应于简并 eval 的 evec 保证是正交的(在任何意义上)。我可以看到如何解决这个问题(在子空间中对角化)，但目前我无法生成这样的示例，而且我太累了，无论如何都无法编写代码。我周末看看是否可以。
关于julia - 如何找到一个复矩阵 A 使得 A * Transpose(A) = B，其中 B = Transpose(B) 并且是复数？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/76884871/