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