Alternating Algorithm
We want to solve the right definite two-parameter eigenvalue
\[\begin{split}A_1u+\lambda B_1u+\mu C_1u=0,\\
A_2v+\lambda B_2v+\mu C_2v=0.\end{split}\]
First we generate matrices satisfying the assumptions, i.e., the matrices \(C_2\) and \(C_1\otimes B_2-C_2\otimes B_1\) are positive definite and \(C_1\) is negative definite.
[1]:
using LinearAlgebra
n=30
A1=randn(n,n); A1=A1'+A1
A2=randn(n,n); A2=A2'+A2
S1=randn(n,n);S2=randn(n,n)
B1=S1*Matrix(Diagonal(rand(n).-0.5))*S1'
B2=S2*Matrix(Diagonal(rand(n).-1.5))*S2'
C1=-S1*S1'
C2=S2*S2';
We can run the alternating algorithm for these matrices.
[2]:
function twopareig_alt(A1,B1,C1,A2,B2,C2,i,j,its)
#input matrices of the two-parameter eigenvalue problem satisfying the definiteness assumptions and desired index
#output eigenvalue (lambda,mu) and eigenvector (u,v) of the desired index
#first we need the dimensions of the matrices
n=size(A1)[1]
m=size(A2)[1]
#we can start with a random vector
u=randn(n)
v=randn(m)
lambda=0
mu=0
for it in 1:its
#compute coefficients for the generalized eigenvalue problem
a=u'*A1*u
b=u'*B1*u
c=u'*C1*u
lhs=a*C2-c*A2
rhs=c*B2-b*C2
l,V=eigen(lhs,rhs)
#asking for the sign of c tells us what eigenvalue to look for to get one of the right index
#if the matrices satsify the assumptions this can be skipped
if c<0
v=V[:,j]
lambda=l[j]
else
v=V[:,m+1-j]
lambda=l[m+1-j]
end
#we can compute mu here
mu=-(a+b*lambda)/c
#the same but for v
a=v'*A2*v
b=v'*B2*v
c=v'*C2*v
lhs=c*A1-a*C1
rhs=b*C1-c*B1
l,U=eigen(lhs,rhs)
if c>0
u=U[:,i]
lambda=l[i]
else
u=U[:,n+1-i]
lambda=l[n+1-i]
end
mu=-(a+b*lambda)/c
end
return lambda,mu,u,v
end
[2]:
twopareig_alt (generic function with 1 method)
We can check this algorithm.
[3]:
for i in 1:n
for j in 1:n
lambda,mu,u,v=twopareig_alt(A1,B1,C1,A2,B2,C2,i,j,5)
#if the eigenvalue has the right index, then the i-/j-th smallest eigenvalue of A+lambda B+mu C is zero
test_i=eigvals(A1+lambda*B1+mu*C1)[i]
test_j=eigvals(A2+lambda*B2+mu*C2)[j]
if abs(test_i)>1e-9 || abs(test_j)>1e-12
println("For index ($i,$j)) we have the eigenvalue errors $(test_i) and $(test_j)")
end
#We can also check the eigenvectors:
u/=norm(u)
v/=norm(v)
test_i=norm((A1+lambda*B1+mu*C1)*u)
test_j=norm((A2+lambda*B2+mu*C2)*v)
if abs(test_i)>1e-9 || abs(test_j)>1e-12
println("For index ($i,$j)) we have the residual norms $(test_i) and $(test_j)")
end
end
end
For index (1,1)) we have the eigenvalue errors -2.299916338630789e-12 and -1.4848391177983952e-12
For index (1,1)) we have the residual norms 1.007565436068525e-10 and 1.446959527213292e-11
For index (1,2)) we have the residual norms 1.7651240588569298e-12 and 2.768298686231474e-12
For index (1,3)) we have the residual norms 1.4927149015939952e-12 and 2.4495324783729644e-12
For index (1,4)) we have the residual norms 1.450088366871532e-12 and 5.047718549398743e-12
For index (1,5)) we have the residual norms 1.409891692982889e-12 and 3.058063692818418e-12
For index (1,6)) we have the residual norms 1.0930788392527906e-12 and 2.415168623843043e-12
For index (1,7)) we have the residual norms 1.0237747393400498e-12 and 6.358159633868717e-12
For index (1,8)) we have the residual norms 1.013481948562707e-12 and 4.1878980054087244e-12
For index (1,9)) we have the residual norms 1.2543936291098397e-12 and 4.876055458818501e-12
For index (1,10)) we have the residual norms 1.1358851444511637e-12 and 5.686467357631265e-12
For index (1,11)) we have the residual norms 9.199832931320417e-13 and 4.60146912621372e-12
For index (1,12)) we have the residual norms 9.911460572124554e-13 and 2.7316725524115034e-12
For index (1,13)) we have the residual norms 1.2144121459219417e-12 and 4.915261707781536e-12
For index (1,14)) we have the residual norms 1.2694459207837628e-12 and 4.557898977387573e-12
For index (1,15)) we have the residual norms 1.2539613973503097e-12 and 3.0961965476899864e-12
For index (1,16)) we have the residual norms 9.075437830562211e-13 and 4.099313738519332e-12
For index (1,17)) we have the residual norms 1.5641695120249347e-12 and 5.12111455850182e-12
For index (1,18)) we have the residual norms 1.1497053033821556e-12 and 2.973476175993261e-12
For index (1,19)) we have the residual norms 1.0457041149040424e-12 and 3.0737406807767206e-12
For index (1,20)) we have the residual norms 1.4965304292341277e-12 and 4.826446567190621e-12
For index (1,21)) we have the residual norms 1.1295721578067586e-12 and 4.4918168353863325e-12
For index (1,22)) we have the residual norms 1.118203667437359e-12 and 3.2179518339844782e-12
For index (1,23)) we have the residual norms 1.4068634023863164e-12 and 3.2819617989152e-12
For index (1,24)) we have the residual norms 1.13316915422272e-12 and 4.2282913491255016e-12
For index (1,25)) we have the residual norms 1.2081529184876352e-12 and 4.2048863964448995e-12
For index (1,26)) we have the residual norms 1.418269445643421e-12 and 2.822022041248379e-12
For index (1,27)) we have the residual norms 1.3739211742128354e-12 and 2.073512774897202e-12
For index (1,28)) we have the residual norms 9.9289988024006e-13 and 3.1499547846239353e-12
For index (1,29)) we have the residual norms 1.0063569313883062e-12 and 3.9122397609116744e-12
For index (1,30)) we have the residual norms 3.767180183681899e-12 and 2.1699173113668634e-11
For index (2,1)) we have the residual norms 3.130155519512396e-11 and 1.3103633329549639e-11
For index (3,1)) we have the residual norms 1.2991383528606218e-11 and 1.1265614181662251e-11
For index (4,1)) we have the eigenvalue errors 5.280326711335341e-13 and -1.8111385211541339e-12
For index (4,1)) we have the residual norms 2.1080591456234276e-11 and 1.0629056734718245e-11
For index (5,1)) we have the residual norms 2.0620636437424063e-11 and 1.498460507732288e-11
For index (6,1)) we have the eigenvalue errors 1.6496797316265218e-12 and -3.1305171330730298e-12
For index (6,1)) we have the residual norms 1.947997650568234e-11 and 1.0286914204963599e-11
For index (7,1)) we have the residual norms 3.100840741967658e-11 and 1.28642403421094e-11
For index (8,1)) we have the residual norms 1.5531744218062812e-11 and 1.2584403918360782e-11
For index (9,1)) we have the eigenvalue errors -2.5236428921575783e-13 and 1.2660243951926684e-12
For index (9,1)) we have the residual norms 2.654641964889775e-11 and 9.936784560451243e-12
For index (10,1)) we have the residual norms 1.7338611349849695e-11 and 9.8552584821935e-12
For index (11,1)) we have the residual norms 1.7860428757969062e-11 and 9.906357160764383e-12
For index (11,30)) we have the residual norms 1.086864138338715e-12 and 1.017237300901148e-12
For index (12,1)) we have the residual norms 2.2495941516704017e-11 and 8.081423331829523e-12
For index (13,1)) we have the residual norms 2.7945089861779047e-11 and 1.735507432054704e-11
For index (14,1)) we have the residual norms 1.454717239910718e-11 and 1.1899834487268948e-11
For index (15,1)) we have the residual norms 1.9431927582512018e-11 and 9.383160063693149e-12
For index (16,1)) we have the eigenvalue errors 2.2618129651975783e-12 and 1.9860648007917634e-12
For index (16,1)) we have the residual norms 2.1156153140884404e-11 and 7.623809002271582e-12
For index (17,1)) we have the eigenvalue errors -8.811332825632847e-14 and 1.1995717477434425e-12
For index (17,1)) we have the residual norms 1.8604859095089048e-11 and 5.456079959925958e-12
For index (18,1)) we have the residual norms 1.9432330591260774e-11 and 9.720040884797248e-12
For index (18,30)) we have the residual norms 1.7737456263244061e-12 and 1.0386888720751964e-12
For index (19,1)) we have the residual norms 1.0372773504751494e-11 and 1.1123797859269707e-11
For index (20,1)) we have the eigenvalue errors -2.1639335812694326e-12 and 1.5648439523404147e-12
For index (20,1)) we have the residual norms 1.9685500998259324e-11 and 7.356739333925452e-12
For index (21,1)) we have the residual norms 1.3831184327542312e-11 and 7.154050408689777e-12
For index (22,1)) we have the eigenvalue errors -2.505918335287778e-13 and -1.8081747430130817e-12
For index (22,1)) we have the residual norms 1.551086352292243e-11 and 9.993258172230044e-12
For index (23,1)) we have the eigenvalue errors 1.4416459696050617e-12 and -1.1581376329330924e-12
For index (23,1)) we have the residual norms 2.731911695355464e-11 and 2.381533651407074e-11
For index (24,1)) we have the eigenvalue errors -2.445060488896789e-12 and 2.4516294371648345e-12
For index (24,1)) we have the residual norms 1.929393593437462e-11 and 7.706650806391231e-12
For index (25,1)) we have the residual norms 1.3564600460091553e-11 and 7.495421196783329e-12
For index (26,1)) we have the eigenvalue errors -2.889263212713911e-12 and -4.7911481208838765e-12
For index (26,1)) we have the residual norms 1.497795528619082e-11 and 1.3095258403830984e-11
For index (27,1)) we have the residual norms 1.0723141379276691e-11 and 1.0106101685834953e-11
For index (28,1)) we have the eigenvalue errors 1.353498144331773e-12 and -3.4596407432218628e-12
For index (28,1)) we have the residual norms 1.505062902930514e-11 and 1.2281943531587308e-11
For index (28,30)) we have the residual norms 1.6470579716265897e-12 and 1.0884893780892688e-12
For index (29,1)) we have the eigenvalue errors 7.552169053861178e-13 and -1.2043763543287786e-12
For index (29,1)) we have the residual norms 9.405523465147644e-12 and 8.577262747148378e-12
For index (29,30)) we have the residual norms 2.2550380008367894e-12 and 1.5189793418324904e-12
For index (30,1)) we have the residual norms 1.5210523409429953e-11 and 6.226523950992213e-11
For index (30,2)) we have the residual norms 3.3850948718604005e-12 and 8.721934557224366e-12
For index (30,3)) we have the residual norms 2.9262738790630724e-12 and 1.1904336834405298e-11
For index (30,4)) we have the residual norms 3.6262209327229847e-12 and 9.928434997868014e-12
For index (30,5)) we have the residual norms 3.3121772620065316e-12 and 6.311290920473303e-12
For index (30,6)) we have the residual norms 3.342618389139544e-12 and 1.09261765387093e-11
For index (30,7)) we have the residual norms 2.41603599118979e-12 and 1.1595960608813501e-11
For index (30,8)) we have the residual norms 2.83160641022739e-12 and 1.535568461606217e-11
For index (30,9)) we have the residual norms 3.4048085423256157e-12 and 1.5966135300251385e-11
For index (30,10)) we have the residual norms 2.733747961179074e-12 and 1.5147696067427455e-11
For index (30,11)) we have the residual norms 2.4482476629546796e-12 and 1.179533665734674e-11
For index (30,12)) we have the residual norms 3.264365875483205e-12 and 1.1972623617852602e-11
For index (30,13)) we have the residual norms 2.999858080296222e-12 and 1.9941247185239848e-11
For index (30,14)) we have the residual norms 3.161991017672155e-12 and 1.2247512189427042e-11
For index (30,15)) we have the residual norms 2.694805338874405e-12 and 1.8458925941025874e-11
For index (30,16)) we have the residual norms 2.719595359557657e-12 and 2.0670152333758558e-11
For index (30,17)) we have the residual norms 3.1422145793905673e-12 and 3.073670781368437e-11
For index (30,18)) we have the residual norms 3.9395990562907766e-12 and 3.1218451815194934e-11
For index (30,19)) we have the residual norms 3.4636867669853755e-12 and 1.830139831007484e-11
For index (30,20)) we have the residual norms 4.066721376410742e-12 and 3.539461080983034e-11
For index (30,21)) we have the residual norms 3.891543382300403e-12 and 2.7970121327029345e-11
For index (30,22)) we have the residual norms 3.53811043435239e-12 and 2.570328724182112e-11
For index (30,23)) we have the residual norms 3.5964712099654897e-12 and 1.7869300209728357e-11
For index (30,24)) we have the residual norms 3.772975236270635e-12 and 1.9313645656251936e-11
For index (30,25)) we have the eigenvalue errors -2.4005950489848034e-13 and 1.1056836896375626e-12
For index (30,25)) we have the residual norms 3.725616175432284e-12 and 2.2184011604106884e-11
For index (30,26)) we have the residual norms 3.530811186056904e-12 and 2.645813381337809e-11
For index (30,27)) we have the residual norms 2.4394909677818813e-12 and 2.7717420586929758e-11
For index (30,28)) we have the eigenvalue errors -4.79398514478226e-13 and 1.285041726083163e-12
For index (30,28)) we have the residual norms 4.845229417054226e-12 and 8.59682055999377e-12
For index (30,29)) we have the residual norms 4.387599954972764e-12 and 8.545348791809817e-12
For index (30,30)) we have the residual norms 4.52599687776896e-12 and 8.280345218977023e-12
These values are all close to machine precision.
If the matrices do not satisfy these assumptions, we can perform affine transformations to the eigenvalue problem. To test this, let us perform some affine transformation to make \(C_1\) and \(C_2\) no longer positive definite.
[4]:
C1=C1+B1; C2=C2+B2;
hcat(eigvals(C1),eigvals(C2))
[4]:
30×2 Array{Float64,2}:
-104.634 -25.3977
-97.3922 -21.675
-83.1846 -17.2486
-75.5028 -13.3026
-63.9637 -11.3747
-59.8023 -10.3087
-53.8328 -7.45019
-51.1422 -6.28786
-48.6185 -6.09173
-38.7105 -3.46527
-33.3855 -2.43435
-31.1964 -1.86886
-26.4929 -1.29252
⋮
-11.5626 0.0436495
-9.62568 0.159252
-9.46295 1.12067
-6.8684 1.27471
-4.5231 2.56337
-3.39391 3.68401
-3.12846 4.59525
-2.21864 6.62235
-0.989113 9.05409
-0.339664 12.058
-0.0757441 19.1556
-0.0101906 25.9177
Now \(C_2\) might not be definite. The algorithm still works:
[5]:
for i in 1:n
for j in 1:n
lambda,mu,u,v=twopareig_alt(A1,B1,C1,A2,B2,C2,i,j,5)
#if the eigenvalue has the right index, then the i-/j-th smallest eigenvalue of A+lambda B+mu C is zero
test_i=eigvals(A1+lambda*B1+mu*C1)[i]
test_j=eigvals(A2+lambda*B2+mu*C2)[j]
if abs(test_i)>1e-9 || abs(test_j)>1e-12
println("For index ($i,$j)) we have the eigenvalue errors $(test_i) and $(test_j)")
end
#We can also check the eigenvectors:
u/=norm(u)
v/=norm(v)
test_i=norm((A1+lambda*B1+mu*C1)*u)
test_j=norm((A2+lambda*B2+mu*C2)*v)
if abs(test_i)>1e-9 || abs(test_j)>1e-12
println("For index ($i,$j)) we have the residual norms $(test_i) and $(test_j)")
end
end
end
For index (1,1)) we have the residual norms 3.355453603877696e-10 and 1.1212841657537528e-11
For index (1,2)) we have the residual norms 2.1343458274495683e-12 and 1.4739619919610336e-12
For index (1,3)) we have the residual norms 1.5160589688392827e-12 and 1.9858340727372666e-12
For index (1,4)) we have the residual norms 1.45367893604648e-12 and 1.6618233596782185e-12
For index (1,5)) we have the residual norms 1.475980712854762e-12 and 1.2727668031506047e-12
For index (1,6)) we have the residual norms 1.0144526667852147e-12 and 1.3794460275459853e-12
For index (1,7)) we have the residual norms 1.6304418281178605e-12 and 1.27321452992426e-12
For index (1,10)) we have the residual norms 1.7053158653190288e-12 and 1.1668220674449482e-12
For index (1,13)) we have the residual norms 8.708985730475105e-12 and 5.872939150745087e-12
For index (1,15)) we have the residual norms 1.4929230562633379e-12 and 2.123988650018873e-12
For index (1,16)) we have the residual norms 1.004247141766847e-12 and 2.8592485600494842e-12
For index (1,17)) we have the residual norms 9.649882772932592e-13 and 4.481109141570219e-12
For index (1,18)) we have the residual norms 1.3079445994386237e-12 and 1.224839099813035e-12
For index (1,19)) we have the residual norms 1.889782850623861e-12 and 1.1035015496199143e-12
For index (1,20)) we have the residual norms 1.0677084036153928e-12 and 4.000306525992516e-12
For index (1,21)) we have the residual norms 1.3401086537464103e-12 and 4.226839612134357e-12
For index (1,22)) we have the residual norms 1.0573002383617575e-12 and 1.4486682844314212e-12
For index (1,23)) we have the residual norms 1.2007167284434399e-12 and 2.1166596667526605e-12
For index (1,24)) we have the residual norms 1.324339325019499e-12 and 4.394704666645571e-12
For index (1,25)) we have the residual norms 1.1894858814343157e-12 and 4.365182992237475e-12
For index (1,26)) we have the residual norms 1.0250065997180808e-12 and 2.520464812330524e-12
For index (1,27)) we have the residual norms 1.3968636910149112e-12 and 2.574822536769229e-12
For index (1,28)) we have the residual norms 1.1665725685650466e-12 and 2.652894769261128e-12
For index (1,29)) we have the residual norms 1.2432766600375633e-12 and 3.688953880111645e-12
For index (1,30)) we have the residual norms 1.79671481653022e-11 and 2.056892662152906e-11
For index (2,1)) we have the eigenvalue errors -1.231887165217819e-11 and -5.2867612840865426e-12
For index (2,1)) we have the residual norms 1.8833268631319488e-10 and 1.3612886230883193e-11
For index (3,1)) we have the residual norms 1.053034644886837e-10 and 9.770569788899569e-12
For index (4,1)) we have the eigenvalue errors -2.70532013223187e-12 and -3.747221063173405e-12
For index (4,1)) we have the residual norms 1.7115209604059094e-10 and 1.1486083211048086e-11
For index (5,1)) we have the eigenvalue errors 8.36171290319404e-12 and -1.9939088361437478e-12
For index (5,1)) we have the residual norms 1.0166527751957823e-10 and 7.524167822641784e-12
For index (6,1)) we have the residual norms 2.714149758036015e-10 and 7.569454818030715e-12
For index (7,1)) we have the eigenvalue errors -6.285036473119984e-12 and -1.5078055872741903e-12
For index (7,1)) we have the residual norms 4.2188259099620515e-10 and 4.9516111788530194e-12
For index (7,25)) we have the residual norms 1.6007909876748411e-10 and 3.5314796338411007e-12
For index (8,1)) we have the residual norms 2.4663905524882437e-10 and 9.039530716299582e-12
For index (9,1)) we have the eigenvalue errors -7.2652637779628146e-12 and -1.0653804181298513e-12
For index (9,1)) we have the residual norms 9.089889240254961e-10 and 1.0196426899699892e-11
For index (9,27)) we have the residual norms 2.0744076614645364e-11 and 2.722273589250362e-12
For index (10,1)) we have the eigenvalue errors 2.497331394197392e-11 and 2.25908347286912e-12
For index (10,1)) we have the residual norms 1.8076773202961229e-10 and 1.0320371519060979e-11
For index (10,27)) we have the residual norms 1.01301586706828e-9 and 3.893030529199728e-12
For index (11,1)) we have the residual norms 3.083139819373278e-10 and 8.554300101869664e-12
For index (11,27)) we have the residual norms 1.2524820065037836e-10 and 3.220965065555092e-12
For index (12,1)) we have the eigenvalue errors -8.820735492594171e-12 and -2.747402628245684e-12
For index (12,1)) we have the residual norms 6.72687753822128e-10 and 8.558540343677372e-12
For index (13,1)) we have the residual norms 7.075842322691325e-10 and 1.042172737776808e-11
For index (14,1)) we have the eigenvalue errors 1.3359514229864314e-11 and -3.010142060931302e-12
For index (14,1)) we have the residual norms 9.470220782887363e-10 and 1.1637694428522675e-11
For index (15,1)) we have the eigenvalue errors -1.0382239382953245e-11 and 1.6463551177670443e-12
For index (15,1)) we have the residual norms 7.402106700568608e-10 and 8.096861289519464e-12
For index (16,1)) we have the eigenvalue errors -2.8391437250115225e-14 and 1.1090520811746811e-12
For index (16,1)) we have the residual norms 4.3113334548188286e-10 and 8.292101462015736e-12
For index (17,1)) we have the eigenvalue errors -6.841393711752888e-13 and 1.2597808441004944e-12
For index (17,1)) we have the residual norms 9.89440436076533e-10 and 8.978047688101634e-12
For index (18,1)) we have the eigenvalue errors -8.422201476163325e-12 and 3.1191936912008806e-12
For index (18,1)) we have the residual norms 8.870092667331252e-10 and 1.5400715071718988e-11
For index (19,1)) we have the residual norms 3.136208629445142e-10 and 1.3253711576702747e-11
For index (20,1)) we have the eigenvalue errors -1.1062423765906313e-11 and 1.0357704899642632e-12
For index (20,1)) we have the residual norms 2.5286509904951474e-10 and 1.0658392091877965e-11
For index (21,1)) we have the eigenvalue errors 1.296407347945506e-11 and -1.1650842452720783e-12
For index (21,1)) we have the residual norms 8.304658501048288e-10 and 7.722034649029588e-12
For index (21,30)) we have the residual norms 1.6538907884440498e-11 and 1.0066521895578263e-12
For index (22,1)) we have the eigenvalue errors -9.546000180667566e-12 and -1.5797331767774984e-12
For index (22,1)) we have the residual norms 4.604965912337181e-10 and 9.679348043067353e-12
For index (23,1)) we have the eigenvalue errors 1.3909160289364049e-11 and -2.0949459333740564e-12
For index (23,1)) we have the residual norms 5.774853285899246e-10 and 9.984493960196829e-12
For index (23,30)) we have the residual norms 1.2692822816226535e-11 and 1.0282405928728515e-12
For index (24,1)) we have the residual norms 5.878666208384777e-10 and 1.456385086886124e-11
For index (25,1)) we have the eigenvalue errors 6.160358068151179e-11 and 1.288910137051153e-12
For index (25,1)) we have the residual norms 4.932237335126846e-10 and 1.8255119615113915e-11
For index (26,1)) we have the eigenvalue errors 1.7870324372374176e-11 and -3.6698648073386485e-12
For index (26,1)) we have the residual norms 6.459594186619013e-10 and 1.0797343914329853e-11
For index (26,30)) we have the residual norms 8.649127148715754e-12 and 1.0341845170464736e-12
For index (27,1)) we have the eigenvalue errors -3.365293886998761e-11 and 1.5018301437704676e-12
For index (27,1)) we have the residual norms 7.900738301510364e-10 and 9.049803118870258e-12
For index (28,1)) we have the residual norms 3.6517926964066404e-10 and 8.469788814277582e-12
For index (28,19)) we have the residual norms 1.5673450151132983e-12 and 1.3025761989925737e-12
For index (28,30)) we have the residual norms 3.5816919709321004e-12 and 1.259747635364954e-12
For index (29,1)) we have the eigenvalue errors -2.4625257381556676e-11 and 1.5459239259884653e-12
For index (29,1)) we have the residual norms 9.282031786196675e-10 and 1.2427184367166689e-11
For index (29,2)) we have the residual norms 2.421085024305741e-11 and 7.903959968097655e-12
For index (29,30)) we have the residual norms 3.765410203173944e-12 and 1.4342670772473637e-12
For index (30,1)) we have the residual norms 1.091880941966486e-9 and 5.661980543492355e-10
For index (30,2)) we have the residual norms 3.493835161321706e-12 and 3.21737125047479e-12
For index (30,3)) we have the residual norms 3.054375206922523e-12 and 3.3079775696351804e-12
For index (30,4)) we have the residual norms 2.7290877223199424e-12 and 3.739465736824038e-12
For index (30,5)) we have the residual norms 2.5048243962893925e-12 and 3.980440655719155e-12
For index (30,6)) we have the residual norms 2.8145117077328435e-12 and 8.140026063708109e-12
For index (30,7)) we have the residual norms 2.4858716133917495e-12 and 3.1663455089733647e-12
For index (30,8)) we have the residual norms 2.6187837909714703e-12 and 4.384670864212289e-12
For index (30,9)) we have the residual norms 2.4126718979098923e-12 and 5.1075313004066325e-12
For index (30,10)) we have the residual norms 3.9046052774716205e-12 and 4.01705572886527e-12
For index (30,11)) we have the residual norms 3.1353981393105137e-12 and 3.721735812604121e-12
For index (30,12)) we have the residual norms 3.4705285065288826e-12 and 3.155262875577935e-12
For index (30,13)) we have the residual norms 2.9025344617461193e-12 and 5.7003241329400656e-12
For index (30,14)) we have the residual norms 3.197525206498398e-12 and 2.6223552847365324e-12
For index (30,15)) we have the residual norms 2.9347866514996045e-12 and 3.903702815995106e-12
For index (30,16)) we have the residual norms 2.9315124943448458e-12 and 6.7075827587362195e-12
For index (30,17)) we have the residual norms 2.8926513520404023e-12 and 7.608073560726328e-12
For index (30,18)) we have the residual norms 2.6951729875872367e-12 and 6.154557874023407e-12
For index (30,19)) we have the residual norms 2.7529198651243984e-12 and 6.549228203095949e-12
For index (30,20)) we have the residual norms 3.940188515456731e-12 and 8.28450367062146e-12
For index (30,21)) we have the residual norms 3.525129909070707e-12 and 6.8402269144736195e-12
For index (30,22)) we have the residual norms 3.599332060020759e-12 and 8.987527338922007e-12
For index (30,23)) we have the residual norms 3.3356349078844303e-12 and 6.8531421104583365e-12
For index (30,24)) we have the residual norms 4.481384641306974e-12 and 7.548479312783976e-12
For index (30,25)) we have the residual norms 3.903732395367044e-12 and 3.611341755224064e-12
For index (30,26)) we have the residual norms 2.0956313370675473e-12 and 6.755710896204468e-12
For index (30,27)) we have the residual norms 3.419578731498809e-12 and 1.2144004777398252e-11
For index (30,28)) we have the residual norms 4.992835717482345e-12 and 5.891036987035551e-12
For index (30,29)) we have the residual norms 3.7686506163673485e-12 and 5.7463910999148395e-12
For index (30,30)) we have the eigenvalue errors 5.030272764894132e-13 and 1.9019417222860556e-12
For index (30,30)) we have the residual norms 4.3272519300009575e-12 and 9.084096441786275e-12
We can perform the steps of Remark 2 to get a problem satisfying the definiteness assumptions. In the same steps we can also diagonolize the matrices \(B_1,B_2,C_1,C_2\). This requires checking multiple cases.
[6]:
function affine_transformation_diagonalization(B1,B2,C1,C2)
#returns diagonal versions of B1,B2,C1,C2 and the corresponding change of basis and affine transformation
l1,U1=eigen(C1)
l2,U2=eigen(C2)
S1=U1;S2=U2
b1=l1;b2=l2;c1=l1;c2=l2
#first diagonilize C1 and C2
#S1,S2 denotes the change of basis
newB1=S1'*B1*S1
newB2=S2'*B2*S2
#we need to remember the affine transformations
T=Matrix(Diagonal(ones(2)))
#check definiteness of C2
if minimum(l2)*maximum(l2)>0
#C2 is defintite
if minimum(l1)*maximum(l1)>0
#C1 is also definite so we only need to diagonalize the B and C matrices
if l1[1]>0
S1=S1*Diagonal(sqrt.(l1.^(-1)))
newB1=Diagonal(sqrt.(l1.^(-1)))*newB1*Diagonal(sqrt.(l1.^(-1)))
b1,U=eigen(newB1)
S1=S1*U
c1.=1
println(norm(S1'*C1*S1-I))
else
S1=S1*Diagonal(sqrt.(-l1.^(-1)))
newB1=Diagonal(sqrt.(-l1.^(-1)))*newB1*Diagonal(sqrt.(-l1.^(-1)))
b1,U=eigen(newB1)
S1=S1*U
c1.=-1
end
if l2[1]>0
S2=S2*Diagonal(sqrt.(l2.^(-1)))
newB2=Diagonal(sqrt.(l2.^(-1)))*newB2*Diagonal(sqrt.(l2.^(-1)))
b2,U=eigen(newB2)
S2=S2*U
c2.=1
else
S2=S2*Diagonal(sqrt.(-l2.^(-1)))
newB2=Diagonal(sqrt.(-l2.^(-1)))*newB2*Diagonal(sqrt.(-l2.^(-1)))
b2,U=eigen(newB2)
S2=S2*U
c2.=-1
end
else
#We will proceed similarly to the second part of the remark
#compute the eigenvalues of the gep B2 u=lambda C2 u
if l2[1]>0
S2=S2*Diagonal(sqrt.(l2.^(-1)))
newB2=Diagonal(sqrt.(l2.^(-1)))*newB2*Diagonal(sqrt.(l2.^(-1)))
l2.=1
else
S2=S2*Diagonal(sqrt.(-l2.^(-1)))
newB2=Diagonal(sqrt.(-l2.^(-1)))*newB2*Diagonal(sqrt.(-l2.^(-1)))
l2.=-1
end
l,U=eigen(newB2)
S2=S2*U
println(norm(S2'*B2*S2-Diagonal(l)))
#We know that C1-1/lambda B1 is definite for any eigenvalue lambda of B2
#C2-1/lambda B2 is semidefinite for extreme eigenvalues of B2
c=10
lambda=l[1]
l3=eigvals(Diagonal(l1)*(lambda-c)-newB1*l2[1])
while minimum(l3)*maximum(l3)<0
c/=2
l3=eigvals(Diagonal(l1)*(lambda-c)-newB1*l2[1])
end
T[2,1]=-l2[1]
T[2,2]=lambda-c
#we performed an affine transformation
b2=l;c2=(lambda-c)*l2-l*l2[1]
newC1=Diagonal(l1)*(lambda-c)-newB1*l2[1]
#we only need to diagonalize B1 and C1
l,U=eigen(newC1)
S1=S1*U
newB1=U'*newB1*U
if l[1]>0
S1=S1*Diagonal(sqrt.(l.^(-1)))
newB1=Diagonal(sqrt.(l.^(-1)))*newB1*Diagonal(sqrt.(l.^(-1)))
l.=1
else
S1=S1*Diagonal(sqrt.(-l.^(-1)))
newB1=Diagonal(sqrt.(-l.^(-1)))*newB1*Diagonal(sqrt.(-l.^(-1)))
l.=-1
end
c1=l
b1,U=eigen(newB1)
S1=S1*U
end
else
#C2 is indefinite, so we can resume as in the remark
#first we can perform an affine transformation to make B1 definite
m=length(l2)
T[1,2]=-newB2[m,m]/l2[m]
newB1+=Diagonal(l1)* T[1,2]
newB2+=Diagonal(l2)* T[1,2]
#we next diagonalize B1
l,U=eigen(newB1)
if l[1]>0
S1=S1*U*Diagonal(sqrt.(l.^(-1)))
newC1=Diagonal(sqrt.(l.^(-1)))*U'*Diagonal(l1)*U*Diagonal(sqrt.(l.^(-1)))
b1.=1
else
S1=S1*U*Diagonal(sqrt.(-l.^(-1)))
newC1=Diagonal(sqrt.(-l.^(-1)))*U'*Diagonal(l1)*U*Diagonal(sqrt.(-l.^(-1)))
b1.=-1
end
#now we diagonalize C1 and perform the necessary affine transformations
l,U=eigen(newC1)
S1=S1*U
lambda=maximum(l)
c=10
l3=eigvals(Diagonal(l2)-newB2*(lambda+c)/b1[1])
while minimum(l3)*maximum(l3)<0
c/=2
l3=eigvals(Diagonal(l2)-newB2*(lambda+c)/b1[1])
end
T[2,:]+=-(lambda+c)/b1[1]*T[1,:]
c1=l-b1*(lambda+c)/b1[1]
l,U=eigen(Diagonal(l2)-newB2*(lambda+c)/b1[1])
S2=S2*U
newB2=U'*newB2*U
if l[1]<0
c2.=-1
S2*=Diagonal(sqrt.(-l.^(-1)))
newB2=Diagonal(sqrt.(-l.^(-1)))*newB2*Diagonal(sqrt.(-l.^(-1)))
else
c2.=1
S2*=Diagonal(sqrt.(l.^(-1)))
newB2=Diagonal(sqrt.(l.^(-1)))*newB2*Diagonal(sqrt.(l.^(-1)))
end
b2,U=eigen(newB2)
S2=S2*U
end
return b1,b2,c1,c2,T,S1,S2
end
[6]:
affine_transformation_diagonalization (generic function with 1 method)
Since we can diagonalize the matrices \(B_1,B_2,C_1,C_2\), we can utilize this in our algorithm:
[7]:
function twopareig_alt_diag(A1,b1,c1,A2,b2,c2,i,j,its)
#input matrices of the two-parameter eigenvalue problem satisfying the definiteness assumptions and desired index
#b1,c1,b2,c2 are the entries of the diagonal matrices B_1,C_1,B_2,C_2
#output eigenvalue (lambda,mu) and eigenvector (u,v) of the desired index
#first we need the dimensions of the matrices
n=size(A1)[1]
m=size(A2)[1]
#we can start with a random vector
u=randn(n)
v=randn(m)
lambda=0
mu=0
for it in 1:its
#compute coefficients for the generalized eigenvalue problem
a=u'*A1*u
b=u'*Diagonal(b1)*u
c=u'*Diagonal(c1)*u
lhs=a*Diagonal(c2)-c*A2
rhs=c*b2-b*c2
#rhs is definite
if minimum(rhs)<0
rhs=Diagonal(sqrt.(-rhs.^(-1)))
lhs=-rhs*lhs*rhs
else
rhs=Diagonal(sqrt.(rhs.^(-1)))
lhs=rhs*lhs*rhs
end
l,V=eigen(Symmetric(lhs))
#asking for the sign of c tells us what eigenvalue to look for to get one of the right index
#if the matrices satsify the assumptions this can be skipped
if c<0
v=rhs*V[:,j]
lambda=l[j]
else
v=rhs*V[:,m+1-j]
lambda=l[m+1-j]
end
#we can compute mu here
mu=-(a+b*lambda)/c
#the same but for v
a=v'*A2*v
b=v'*Diagonal(b2)*v
c=v'*Diagonal(c2)*v
lhs=c*A1-a*Diagonal(c1)
rhs=b*c1-c*b1
if minimum(rhs)<0
rhs=Diagonal(sqrt.(-rhs.^(-1)))
lhs=-rhs*lhs*rhs
else
rhs=Diagonal(sqrt.(rhs.^(-1)))
lhs=rhs*lhs*rhs
end
l,U=eigen(Symmetric(lhs))
if c>0
u=rhs*U[:,i]
lambda=l[i]
else
u=rhs*U[:,n+1-i]
lambda=l[n+1-i]
end
mu=-(a+b*lambda)/c
end
return lambda,mu,u,v
end
[7]:
twopareig_alt_diag (generic function with 1 method)
We can now effectively compute all eigenvalues of the two-parameter eigenvalue problem.
[8]:
function twopareig_alt(A1,B1,C1,A2,B2,C2,its)
#first we perform affine transformations and diagonalization of the matrices B1,B2,C1,C2
b1,b2,c1,c2,T,S1,S2=affine_transformation_diagonalization(B1,B2,C1,C2)
A1new=Symmetric(S1'*A1*S1)
A2new=Symmetric(S2'*A2*S2)
n=length(b1)
m=length(b2)
lambdas=zeros(n,m)
mus=zeros(n,m)
U=zeros(n,m,n)
V=zeros(n,m,m)
for i in 1:n
for j in 1:m
#compute the eigenvalues of the transformed eigenvalue problem
lambda,mu,u,v=twopareig_alt_diag(A1new,b1,c1,A2new,b2,c2,i,j,its)
#lambda,mu,u,v=twopareig_alt(A1new,Matrix(Diagonal(b1)),Matrix(Diagonal(c1)),A2new,Matrix(Diagonal(b2)),Matrix(Diagonal(c2)),i,j,its)
#we can get the eigenvalues and eigenvectors of the original problem easily
lambdas[i,j]=T[1,1]*lambda+T[2,1]*mu
mus[i,j]=T[1,2]*lambda+T[2,2]*mu
U[i,j,:]=S1*u
V[i,j,:]=S2*v
end
end
lambdas,mus,U,V
end
#we can compare this method to the one without performing affine transformations
function twopareig_alt_compare(A1,B1,C1,A2,B2,C2,its)
n=size(B1)[1]
m=size(B2)[1]
lambdas=zeros(n,m)
mus=zeros(n,m)
U=zeros(n,m,n)
V=zeros(n,m,m)
for i in 1:n
for j in 1:m
#compute the eigenvalues of the transformed eigenvalue problem
lambdas[i,j],mus[i,j],U[i,j,:],V[i,j,:]=twopareig_alt(A1,B1,C1,A2,B2,C2,i,j,its)
end
end
lambdas,mus,U,V
end
[8]:
twopareig_alt_compare (generic function with 1 method)
Diagonalizing the matrices \(B_1,B_2,C_1,C_2\) can save significant time when computing multiple eigenvalues for larger matrices:
[9]:
println()
println("Time of algorithm with diagonalization")
@time lambdas,mus,U,V= twopareig_alt(A1,B1,C1,A2,B2,C2,5)
@time lambdas,mus,U,V= twopareig_alt(A1,B1,C1,A2,B2,C2,5)
println()
finalerror=0
finalresnorms=0
for i in 1:n
for j in 1:n
lambda=lambdas[i,j]
mu=mus[i,j]
test_i=minimum(abs.(eigvals(A1+lambda*B1+mu*C1)))
test_j=minimum(abs.(eigvals(A2+lambda*B2+mu*C2)))
if abs(test_i)>1e-9 || abs(test_j)>1e-9
println("For index ($i,$j)) we have the eigenvalue errors $(test_i) and $(test_j)")
end
finalerror+=abs(test_i)+abs(test_j)
u=U[i,j,:]
v=V[i,j,:]
u/=norm(u)
v/=norm(v)
test_i=norm((A1+lambda*B1+mu*C1)*u)
test_j=norm((A2+lambda*B2+mu*C2)*v)
if abs(test_i)>1e-9 || abs(test_j)>1e-9
println("For index ($i,$j)) we have the residual norms $(test_i) and $(test_j)")
end
finalresnorms+=abs(test_i)+abs(test_j)
end
end
println()
println("the sum of errors is $finalerror")
println("the sum of residual norms is $finalresnorms")
println()
println("Time of algorithm without diagonalization and affine transformation")
@time lambdas,mus,U,V= twopareig_alt_compare(A1,B1,C1,A2,B2,C2,5)
@time lambdas,mus,U,V= twopareig_alt_compare(A1,B1,C1,A2,B2,C2,5)
println()
finalerror=0
finalresnorms=0
for i in 1:n
for j in 1:n
lambda=lambdas[i,j]
mu=mus[i,j]
test_i=minimum(abs.(eigvals(A1+lambda*B1+mu*C1)))
test_j=minimum(abs.(eigvals(A2+lambda*B2+mu*C2)))
if abs(test_i)>1e-9 || abs(test_j)>1e-9
println("For index ($i,$j)) we have the eigenvalue compare errors $(test_i) and $(test_j)")
end
finalerror+=abs(test_i)+abs(test_j)
u=U[i,j,:]
v=V[i,j,:]
u/=norm(u)
v/=norm(v)
test_i=norm((A1+lambda*B1+mu*C1)*u)
test_j=norm((A2+lambda*B2+mu*C2)*v)
if abs(test_i)>1e-9 || abs(test_j)>1e-9
println("For index ($i,$j)) we have the compare residual norms $(test_i) and $(test_j)")
end
finalresnorms+=abs(test_i)+abs(test_j)
end
end
println()
println("the sum of compare errors is $finalerror")
println("the sum of compare residual norms is $finalresnorms")
println()
#we can also compare this with the time we need to solve the assiociated n^2 times n^2 eigenvalue problem
println("Time of computing the associated n^2 times n^2 eigenvalue problem")
@time eigen(kron(A1,C2)-kron(C1,A2),kron(C1,B2)-kron(B1,C2));
@time eigen(kron(A1,C2)-kron(C1,A2),kron(C1,B2)-kron(B1,C2));
Time of algorithm with diagonalization
12.543337 seconds (23.34 M allocations: 1.633 GiB, 5.92% gc time)
2.683528 seconds (263.29 k allocations: 571.902 MiB, 1.39% gc time)
the sum of errors is 1.1980345955734464e-9
the sum of residual norms is 2.161697666562168e-8
Time of algorithm without diagonalization and affine transformation
4.252455 seconds (449.04 k allocations: 690.123 MiB, 1.44% gc time)
3.872379 seconds (393.31 k allocations: 687.333 MiB, 1.48% gc time)
For index (10,27)) we have the compare residual norms 1.1863507944073396e-9 and 9.981404090838363e-12
For index (14,1)) we have the compare residual norms 1.0073722157766232e-9 and 1.02716614281728e-11
For index (17,1)) we have the compare residual norms 1.0755401831257006e-9 and 8.78664527667746e-12
For index (30,1)) we have the compare residual norms 1.1265770822014743e-9 and 4.163646323547098e-11
the sum of compare errors is 8.833694326769068e-10
the sum of compare residual norms is 1.939129171690429e-8
Time of computing the associated n^2 times n^2 eigenvalue problem
5.366581 seconds (32.25 k allocations: 57.604 MiB, 0.10% gc time)
5.071025 seconds (31 allocations: 55.936 MiB, 0.10% gc time)
Diagonalizing did not increase the error! But for large \(n\) we can save a significant amount of time!