連立一次方程式反復法

参照資料: Chapter 3: Sections 1 （スライド）

1 ヤコビ解法（）

1.1 ヤコビ解法

入力： 行列A,   ベクトルb,  初期ベクトルx0,  相次ぐ二つの近似値の制度prec,  最大反復の回数iter
出力：　Ax=bを満たすxの近似値.

jacobi[A_,b_,x0_,prec_,iter_]:=(Module[{j,y,i,k,norm,n,x},
n=Length[b];
x=x0;
y=Table[0,{i,1,n}];
For[k=1,k≤iter,k++,
norm=0;
For[i=1,i≤n,i++,
y[[i]]=b[[i]];
For[j=1,j≤n,j++,
If[j≠i,
y[[i]]-=A[[i,j]]*x[[j]];
]
];
y[[i]]*=1./A[[i,i]];
If[Abs[y[[i]]-x[[i]]]>norm,norm=Abs[y[[i]]-x[[i]]]];
];
Print[y," , norme:",norm," , k:",k];
If[norm<prec,Return[y]];
x=y;
];
Return["Iterations not enough"];
];
)

A={{10,-1,2,0},{-1,11,-1,3},{2,-1,10,-1},{0,3,-1,8}}
b={6,25,-11,15}

1.2 ガウス・サイデルによる改善。

GaussSeidel[A_,b_,x0_,prec_,iter_]:=(Module[{j,y,i,k,norm,n,x},
n=Length[b];
x=x0;
y=Table[0,{i,1,n}];
For[k=1,k≤iter,k++,
norm=0;
For[i=1,i≤n,i++,
y[[i]]=b[[i]];
For[j=1,j<i,j++,
y[[i]]-=A[[i,j]]*y[[j]];
];
For[j=i+1,j≤n,j++,
y[[i]]-=A[[i,j]]*x[[j]];
];
y[[i]]*=1./A[[i,i]];
If[Abs[y[[i]]-x[[i]]]>norm,norm=Abs[y[[i]]-x[[i]]]];
];
Print[y," , norme:",norm," , k:",k];
If[norm<prec,Return[y]];
x=y;
];
Return["Iterations not enough"];
];
)

上記の行列Aで、jabobi解法とGaussSeidel解法を比べる。

jacobi[A,b,{0.,0.,0.,0.},0.00001,20]

GaussSeidel[A,b,{0.,0.,0.,0.},0.000000001,20]

LinearSolve[A,b]

Another Example (Exercise 7.3.1 in Burden Faires p. 449)

A={{4,1,-1},{-1,3,1},{2,2,5}};c={5,-4,1}

GaussSeidel[B,d,{0,0,0},10^(-4),3]

DA=DiagonalMatrix[{4,3,5}]

UA={{0,-1,1},{0,0,-1},{0,0,0}}

LA=DA-UA-A

Tg=Inverse[DA-LA].UA

CharacteristicPolynomial[Tg,λ]

Eigenvalues[Tg]//N

Tj=Inverse[DA].(LA+UA)

Tj//MatrixForm

CharacteristicPolynomial[Tj,λ]

Eigenvalues[Tj]//N

GaussSeidel[A,c,{0,0,0},10^(-5),5]

jacobi[A,c,{0,0,0},10^(-5),5]

Another Example

B={{-2,1,1/2},{1,-2,-1/2},{0,1,2}}

d={4,-4,0}

jacobi[B,d,{0,0,0},10^(-4),5]

1.3 Convergnece: counter example 収束：反例

B={{1,2,-2},{1,1,1},{2,2,1}}
bb={7,2,5}

Solve B.x=bb

Det[B]

d=DiagonalMatrix[{1,1,1}]
u={{0,-2,2},{0,0,-1},{0,0,0}}
l={{0,0,0},{-1,0,0},{-2,-2,0}}

d-u-l==B

Tj=Inverse[d].(l+u)

Eigenvalues[Tj]

Spectral radius of Tj is zero. Try Gauss-Seidel

Tg=Inverse[d-l].u

Eigenvalues[Tg]

Spectral radius larger than 2, so no convergence. Let’s check

GaussSeidel[B,bb,{0,0,0},0.00001,10]

2 緩和法 Relaxation method

SOR[A_,b_,x0_,prec_,iter_,omega_]:=(Module[{j,y,i,k,norm,n,x},
n=Length[b];
x=x0;
y=Table[0,{i,1,n}];
For[k=1,k≤iter,k++,  (* main iteration *) 
norm=0;
For[i=1,i≤n,i++,
y[[i]]=b[[i]];
For[j=1,j<i,j++,
y[[i]]-=A[[i,j]]*y[[j]];
];
For[j=i+1,j≤n,j++,
y[[i]]-=A[[i,j]]*x[[j]];
];
y[[i]]*=(0.+omega/A[[i,i]]);
y[[i]]+=(1-omega)*x[[i]];
If[Abs[y[[i]]-x[[i]]]>norm,norm=Abs[y[[i]]-x[[i]]]];
];
Print[y," , norme:",norm," , k:",k];
If[norm<prec,Return[y]];
x=y;
];
Return["Iterations not enough"];
];)

A={{4,3,0},{3,4,-1},{0,-1,4}}
b={24,30,-24}

A//MatrixForm

x0={1,1,1}

$MinPrecision=10

SOR[A,b,x0,10^(-6),20,1.25]

SOR[A,b,x0,10^(-7),40,1]

d=DiagonalMatrix[{4,4,4}]

l={{0,0,0},{-3,0,0},{0,1,0}}

u={{0,-3,0},{0,0,1},{0,0,0}}

d-u-l

Map[MatrixForm,{u,l,d,d-l}]

Inverse[d-l]

9/16//N

Inverse[d-5/4* l].{-24,30,-24}

d-(5/4)*l//MatrixForm

-1/4*d+5/4*u//MatrixForm

3 共役勾配法 Conjugate gradient

3.1 Steepest (=Gradient) descent 最速降下　（＝勾配降下）

steepest[A_,b_,C_,x0_,iter_,prec_]:=(Module[{x,r,x1,q,t,numt,ar,denomt,i,k,normr},
x=x0;
k=1;
While[k≤iter,
r=b-A.x//N;
normr=Max[Abs[r]];
If[normr<prec,Print["nb of iteration ",k];Return[x]];
numt=Sum[i^2,{i,r}]; (* numerator t =∥r∥^2 *) 
ar=A.r;
denomt=r.ar; (* denominator t = ∥r((Subscript[∥, A])^2) *)
t=numt/denomt;
x+=t*r;If[k>iter-5,Print["iteration ",k," ", x," norm r=",normr]];
k+=1;
];
Print[x];Return["Not enough iterations"];
])

A={{4,3,0},{3,4,-1},{0,-1,4}}
b={24,30,-24}

steepest[A,b,0,{0,0,0},30,0.001]

slowConjGrad[A_,b_,C_,x0_,iter_,prec_]:=(Module[{x,r,q,t,numt,av,denomt,i,k,normr,vnew,v,gv,n},
x=x0;
k=1;
v={};n=Length[x0];
While[k≤iter,
r=b-A.x//N;
vnew=gramSchmidtAux[r,v,A,False];
v=Append[v,vnew];
normr=Max[Abs[r]];
If[normr<prec,Print["nb of iteration ",k];Return[x]];
numt=vnew.r; (* numerator t =∥r∥^2 *) 
av=A.vnew;
denomt=vnew.av; (* denominator t = ∥r((Subscript[∥, A])^2) *)
t=numt/denomt;
x+=t*vnew;(* If[k>iter-5,Print["iteration ",k," ", Table[x[[i]],{i,{1,Floor[n/4],Floor[n/3],Floor[n/2],Floor[2n/3],n}}]," norm r=",normr]]; *)
Print["iteration ",k," ", x," norm r=",normr, "  v=",vnew, "  t=",t];
k+=1;
];
Return["Not enough iterations"];
])

conjugateGrad[A,b,Cond,x0,3,0.0001]

slowConjGrad[A,b,0,{0,0,0},3,0.0001];

3.2 Conjugate gradient descent 共役勾配（降下）

INPUT: invertible matrix A, vector b
	preconditioned matrix C
	initial approximation x0
	max nbr of iterations iter
	tolerance (precision) prec
OUTPUT: approximation x of the system  Ax=b 
	residual 残差 or1,...r,1

conjugateGrad[A_,b_,C_,x0_,iter_,prec_]:=(Module[{x,r,w,v,z,n,beta,k,i,norme,a,u,t,s},
x=x0;
r=b-A.x//N; (* r0  *)
w=C.r; (* w0 *)
v=C.w; (* v1 *)  
a=Sum[i^2,{i,w}]; (* (norme r)^2*) 
n=Length[x0];
k=1;
norme=0;
For[i=1,i≤n,i++,
If[Abs[v[[i]]]>norme,norme=Abs[v[[i]]]];
];
While[k≤iter,
If[norme<prec,
Return[{x,r}];
];
u=A.v; (* u = A vk*)
t=(a/v.u)//N;  (* tk *)
x=x+t*v;
(* If[k>iter-5,Print["iteration ",k, " ",Table[x[[i]],{i,{1,Floor[n/4],Floor[n/3],Floor[n/2],Floor[2n/3],n}}]]]; (* ," norme r=",norme]; *)  *)
Print["iteration ",k, " ",Table[x[[i]],{i,{1,n}}]]; (* ," norme r=",norme]; *)
r-=t*u; (* r=r-t Av *)
w=C.r;
beta=Sum[i^2,{i,w}];
If[Abs[beta]<prec,
norme=0;
For[i=1,i≤n,i++,If[Abs[r[[i]]]>norme,norme=Abs[r[[i]]]]];
If[norme<prec,Return[{x,r}]];
];
s=beta/a;  (* s = sk *)
v=Transpose[C].w+ s*v; 
a=beta;
k++;
];
If[k>iter,Return["Maximal number of iterations exceeded"]];
])

Exercise 1 of Homework January 5th

A={{4,3},{3,4}}
b={2,5}

conjugateGrad[A,b,IdentityMatrix[2],{0,0},10,0.00001]

slowConjGrad[A,b,0,{0,0},2,0.00001]

steepest[A,b,IdentityMatrix[2],{0,0},10,0.0001]

x1=29/176 * b //N

t1=29/176//N

b-A.x1

r1=b-t1*{23,26}

r1.(A.r1)

t2=Out[48]/Out[49]

x2=x1+t2*r1

Gram-Schmidt

V1=b

r1

{ r1.A.b  , r1.{23,26}}

V2=r1+(22.5511/176)*V1

22.5511/176

V2.A.V1

(* T2 = < V2,r1>/<V2,V2(Subscript[>, A]) *) {num=V2.r1 , den=V2.A.V2 ,T2 =num/den}

x1+T2*V2

With conjugate gradient formula

r1.r1

s1=Out[99]/29

{r0=b,x0={0,0},v1=b,b.b,b.(A.b),t1=b.b/(b.A.b)}//N

{x1=t1*b,r1=r0-t1*A.v1}//N

x0={0,0,0}
Cond=IdentityMatrix[3]

conjugateGrad[A,b,Cond,x0,3,0.0001]

Example page 490 of Burden-Faires. B below is poorly conditioned

B={{0.2,0.1,1,1,0},{0.1,4,-1,1,-1},{1,-1,60,0,-2},{1,1,0,8,4},{0,-1,-2,4,700}}

B//MatrixForm

Eigenvalues[B]  (* poorly conditioned *)

d={1,2,3,4,5}

Poor convergence before last step (fifth step) is NG

conjugateGrad[B,d,IdentityMatrix[5],{0,0,0,0,0},10,0.0001]

M=RandomReal[{-9.,9.},{100,100}];Q=Transpose[M].M;

Eigenvalues[Q]

conjugateGrad[Q,Table[i/(10. +i),{i,1,100}],IdentityMatrix[100],Table[0.0,100],200,0.001]

steepest[Q,Table[i/(10. +i),{i,1,100}],IdentityMatrix[100],Table[0.0,100],100,0.1]

LinearSolve[M,Table[i/(10.+i) , {i,100}]]

Table[1,10]

3.3 Preconditioning 前処理

Test1: Small condition number implies good convergences:
We build positive definite symmetric matrices with small condition number.

gramSchmidtAux[vect_,listVect_,A_,normal_]:=(Module[{iV,n,vect0},
(* orthogonlize vext wrt to the listVect. normal is a boolean -> orthonormalize, or not *)
n=Length[listVect];
If[n==0 ,
If[normal,
Return[vect/N[Sqrt[vect.(A.vect)]]],
Return[vect];
];
];
If[normal,
vect0=vect-Sum[N[vect.(A.iV)] * iV,{iV,listVect}];
Return[vect0/N[Sqrt[vect0.(A.vect0)]]];
,Return[vect-Sum[(N[vect.(A.iV)] / N[iV.(A.iV)])* iV,{iV,listVect}]];
];
];)

gramSchmidt[vectors_,A_,normal_]:=(Module[{n,i,ortho,cur,curVect},
n=Length[vectors];(* Print["n=",n]; *)
ortho={};
For[i=1,i≤n,
cur=vectors[[i]];
(* Print["i=",i," ortho=",ortho," current vector=", cur]; *)
curVect=gramSchmidtAux[cur,ortho,A,normal];
ortho=Append[ortho,curVect];
i++;
];
Return[ortho];
];
)

vect={{1,23423,2121,2342,121},{42342.34,6243,513431,51,1},{53,-3,51,51,16},{1,-2,1,-3,7}}

gramSchmidt[vect,IdentityMatrix[5],True]

Generating positive definite matrices with a small condition number

generatePosDefMat[size_,eig_,modulus_]:=(Module[{D,Q,Q1},
D=DiagonalMatrix[eig];
Q1=Table[Table[RandomReal[{-modulus,modulus}],{i,1,size}],{j,1,size}];
Q=gramSchmidt[Q1,IdentityMatrix[size],True];
Return[Transpose[Q].D.Q];
]);

randomPosDef[size_,min_,max_]:=(generatePosDefMat[size,Union[Table[RandomReal[{min,10*min}],size/10],Table[RandomReal[{10*min+1,max*(0.9)}],size*4/5],Table[RandomReal[{(max*0.9)+1,max}],size/10]],max]);

Union[{1,3,4},{10,11,12},{20,21,22}]

randomPosDef[4,1,3]//MatrixForm

test=randomPosDef[300,0.001,10000];

Eigenvalues[test]

result=conjugateGrad[test,Table[i/(10. +i),{i,1,300}],IdentityMatrix[300],Table[0.0,300],399,0.0001];

Table[result[[1,i]],{i,{1,Floor[n/4],Floor[n/3],Floor[n/2],Floor[2n/3],n}}]

n=300;testsol=LinearSolve[test,Table[i/(10. +i),{i,1,300}]];Table[testsol[[i]],{i,{1,Floor[n/4],Floor[n/3],Floor[n/2],Floor[2n/3],n}}]

resultslow=slowConjGrad[test,Table[i/(10. +i),{i,1,300}],IdentityMatrix[300],Table[0.0,300],300,0.0001];

Table[resultslow[[i]],{i,{1,Floor[n/4],Floor[n/3],Floor[n/2],Floor[2n/3],n}}]