Table of Contents
1 W From the Eigen Space \(Q^{T}\)
Dimension of matrices is \(500\times500\). The next figures show the backward error evolution at each iterations of the CG (blue line, \(k\) = 0) and CG-DEF.
Top labels are the strategies used to set \(k\) the number of eigenvectors and \([S1|S2|S3]\) which part of the spectrum is deflated. Right labels \([1E2|1E4|1E6]\)
are vlues of \(\theta\) the distance between LDC or MDC from centered eigen cluster.
1.1 Matrices without clustered eigenpairs:
When matrices don't have LDC or MDC, deflated only a few eigenvectors from LD don't reduced significantly the number of iterations to reach stopping criterion (fig:AOWG , fig:A1WG). When the deflation use the complete eigenspace, the method converge at the iteration \(0\). In this case, deflation effect appear more like preconditioner apply to linear operator \(A0\) and \(A1\).
Figure 1: Backward error when deflating with \(k\) eigenvectors, \(A0\).
Figure 2: Backward error when deflating with \(k\) eigenvectors, \(A1\).
1.2 Matrices with clustered eigenpairs:
1.2.1 Cluster of size 1:
In case of LDC like in fig:A2WG (S3-LD), deflated the LD eigenvector is enough to reduce the number of iterations. Afterwhat (\(k>1\)) deflation
have the same effect like for matrices without clustered eigenvalues.
If the cluster is close to MD like in fig:A3GW, deflation need to be perform with MD eigenvectors (S3-MD) to speed up the convergence
and if matrix has to cluster (LDC and MDC), deflation can be set up with LDC and MDC but best speed up is obtain when deflating with both LDC and MDC.
1.2.2 Cluster of size 4:
Figure 3: Backward error when deflating with \(k\) eigenvectors, \(A2\).
Figure 4: Backward error when deflating with \(k\) eigenvectors, \(A3\).
Figure 5: Backward error when deflating with \(k\) eigenvectors, \(A4\).
Even with a larger cluster, figA2WG4, figA3WG4, figA4WG4 show that deflatind with eigenvectors from cluster speed up the convergence. When all clustered eigenvectors are used the convergence is linear.
To resume this experimentation, deflation need to used eigenvectors from cluster eigenvectors to reduce number of iterations. We show also that
\(\theta\) do not influence the deflation directly but the condition number. It's why CG-DEF fit well for ill-conditionned matrices.
# Testing CG-DEF with W given : from DDMPY import * import os import sys print("Testing CG-DEF with W=Q.T:") # Init Test cases : mat = ["A0", "A1", "A2", "A3", "A4"] theta = [1e2, 1e4, 1e6] maxrun = 3 part = ["LD", "MD", "LMD"] SI = ['S0', 'S1', 'S2', 'S3'] n = 500 nk = 1 if(len(sys.argv) == 1):n = 500 if(len(sys.argv) == 2):n = int(sys.argv[1]) if(len(sys.argv) == 3): n = int(sys.argv[1]) nk = int(sys.argv[2]) rep = "n_"+ str(n) + "_k_" + str(nk) + "_" + "EXP1_1" Wgiven = True for Imat in mat: # Build Matrix and EigenPairs : wtheta = theta if Imat == "A0" or Imat == "A1": wtheta = [1e2] for Itheta in wtheta: if Itheta == 1e2: Stheta="1E2" if Itheta == 1e4: Stheta="1E4" if Itheta == 1e6: Stheta="1E6" A,DA,QTA = setmatrix4exp(which=Imat, n=n, k=nk, theta=Itheta) for ISI in SI: if ISI == "S3": wpart = part if ISI != "S3": wpart = ["LD"] for Ipart in wpart: if Ipart == "LMD": #nk1 = 2; nk2 = 6; step = 2 IkListe = [2, 4] else: #nk1 = 1; nk2 = 5; step = 1 IkListe = [1, 2, 3, 4] if ISI == 'S0' or ISI == 'S1': #nk1 = 1; nk2 = 2; step = 1 IkListe = [1] if ISI == 'S2' : #nk1 = 1; nk2 = 10; step = 1 IkListe = [1, 2, 3, 4, 10, 50, n] irr = "1"; enum = 0; enum_c = 0 for Ik in IkListe:#range(nk1, nk2, step): k=Ik; l1=0; l2=10; Si=ISI; which=Ipart; RHR="Wgiven"; sigma=1e-7 for Irun in range(1, maxrun+1): # Build exact solutions and right-hand sides : xexa = np.random.rand(A.shape[0],2) b = A@xexa ns = 1 W = None strat = [Wgiven, k, l1, l2, Si, which, ns, RHR, sigma, True] # Solve: for i in range(b.shape[1]): if ns > 1 : enum = enum+1; enum_c = enum if ns == 1 : enum_c = 0 Tags = [Imat, str(Irun), Stheta, rep, irr, str(enum_c), ISI, Ipart] # Init W (deflation Subspace) with EigenSpace : if i>0 and Si != 'S0': if Si == 'S1' : W = QTA[...,0:n].copy() if Si == 'S2' : W = QTA[...,0:k].copy() if Si == 'S3' : if Ipart == 'LD': W = QTA[...,0:k].copy() if Ipart == 'MD': W = QTA[...,n-k:n].copy() if Ipart == 'LMD': W = np.hstack((QTA[...,0:int(k/2)], QTA[...,n-int(k/2):n])) s = ConjGrad(A, ritz=True, tol=1e-16, maxiter=n, true_res=False, xexa=xexa[...,i:i+1], W=W, strat=strat, tags=Tags) irr = "2" if Si != 'S0': x, W = s.dot(b[...,i:i+1]) csv_eigen(A, W, DA, QTA, strat, tags=Tags) else: x = s.dot(b[...,i:i+1]) ns += 1 strat = [Wgiven, k, l1, l2, Si, which, ns, RHR, sigma, True] print("|------> system solved : i=",i) print(s) print("||x-xexa||_{A}=",np.sqrt( (x-xexa[...,i:i+1]).T@A@(x-xexa[...,i:i+1]) ) ) print("||b-A@x||_{2}/||b||_{2}",la.norm(b[...,i:i+1]-A@x)/la.norm(b[...,i:i+1]) ) print("Condition number:") ritz = s.ritz_values() print(np.amax(ritz)/np.amin(ritz)) cmd = "R -q --vanilla < EXP1_1.R" + " --args " + rep + " " + RHR os.system(cmd) cmd = "R -q --vanilla < REXP1_1.R" + " --args " + rep os.system(cmd)