LU-, Cholesky felbontások

Bevezetés

A mátrix alábbi háromszög felbontásait L U illetve (R^H)R felbontásoknak nevezzük, ha
A=L U ((Subscript[l, jj])=1), ahol L alsó, U felső trianguláris;
A=(R^H)R, ahol R felső trianguláris és Ha konjugált transzponált operátor.

MYLUDEMO: Interaktív demonstráció az LU felbontáshoz

Kiindulás: első mátrix: elimináló mátrix (L1), második mátrix: A (rendszermátrix)

Black-Box: LUDecomposition[mat]

MatrixForm[LUDecomposition[(({{1, 0, 9}, {2, 4, 6}, {8, 9, 9}}))][[1]]]

Vegyük észre hogy ebben a mátrixban benne van az LU felbontásra vonatkozó minden információ!

LUDEMO programkód

AUXs

Aux Codes

CheckLUMinors0[A_]:=And@@Table[Minors[A,i][[1,1]]≠0,{i,Dimensions[A][[1]]}]

CheckLUMinors[A_]:=And@@Table[Det[Take[A,i,i]]≠0,{i,Dimensions[A][[1]]}]

Clear[CreateLUMatrix]

CreateLUMatrix[n_]:=Module[{l=True,T},While[l,T=Table[Random[Integer,{0,10}],{n},{n}];If[CheckLUMinors[T]∧LUDecomposition[T][[2]]==Range[n],l=False]];T]

Clear[TM]

TM[A_,n_]:=Module[{II=IdentityMatrix[Dimensions[A][[1]]]},ReplacePart[II,Table[Rule[{j,n},-A[[j,n]]/A[[n,n]]],{j,n+1,Dimensions[A][[1]]}]]];

Clear[TSteps]

TSteps[ML_,n_]/;n==S:=ML;

TSteps[ML_,n_]:=Module[{A},If[n>S,S++;A=Dot@@Last[ML];Append[ML,{TM[A,n],A}],S--;Drop[ML,-1]]]

SM[ML_]:={MatrixForm[ML[[1]]],MatrixForm[ML[[2]]]};

Init[n_]:=Module[{M=CreateLUMatrix[n]},{{TM[M,1],M}}]

S=1;

Main

Main Manipulate

Manipulate[Row[SM[Last[MatrixList=TSteps[MatrixList,s]]]],Style["LUDEMO—2009",Red,16],Button[Init,s=1;S=1;MatrixList=Init[n]],{{n,3},2,4,1},{s,1,n,1},Initialization:>(MatrixList=Init[3]),SaveDefinitions->True]

A) Hasonlítsuk össze A Demo U mátrixát az  LUDecomposition[[1]]//MatrixForm eredményével.
B)  Hasonlítsuk össze  (Ln...L1)^(-1).U mátrixot  A-val
C) Magyarázzuk el a Demo működését

LUDecomposition[(({{1, 0, 7}, {1, 4, 5}, {3, 9, 0}}))][[1]]//MatrixForm

LUDecomposition[(({{3, 1, 1}, {4, 8, 2}, {10, 10, 1}}))][[1]]//MatrixForm

LU (LR) Parketta (implementáció)

(({{■, ■, ■, ■}, {1, ■, ■, ■}, {↓ , ■, ■, ■}, {■, ■, □, □}}))(({{1, ->, ■, ■}, {■, ■, ■, ■}, {■, ■, ■, ■}, {■, ■, □, □}}))

Cél: LU (LR) felbontás, ahol L alsó trianguláris, (Subscript[l, ii])=1
U felülről trianguláris

Megjegyzés . Nem minden mátrixnak létezik ilyen felbontása, még akkor sem, ha az nemszinguláris.
Ha létezik, nem feltétlen egyértelmű.

Black Box: Mathematica LUDecomposition[matrix]
White Box: LU0

Vigyázat! Az egyszerűség kedvéért L, U globális változók.

LU felbontás programkódja

Clear[A];Clear[LU0];

Init

LU0[A_]:=(L=IdentityMatrix[Length[A]];U=DiagonalMatrix[Table[0,{i,Length[A]}]];LU0[A,1])

(Subscript[U, j k]) megh.: (Subscript[a, j k])=(∑_(l=1)^(j-1))(Subscript[L, j l])(Subscript[U, l k])+(Subscript[U, j k])
(Subscript[L , k j]) has.

LU0[A_,j_]/;j>Length[A]:={L,U};
LU0[A_,j_]:=(For[k=j,k≤Length[A],k++,U[[j,k]]=A[[j,k]]-Sum[L[[j,l]] U[[l,k]],{l,j-1}]];If[j<Length[A],For[k=j+1,k≤Length[A],k++,L[[k,j]]=If[U[[j,j]]≠0,(A[[k,j]]-Sum[L[[k,l]] U[[l,j]],{l,j-1}])/U[[j,j]],If[A[[k,j]]-Sum[L[[k,l]] U[[l,j]],{l,j-1}]==0,1,"x"]]]];LU0[A,j+1])

Progamozási Technikák

If
If[cond,Truebranch,Falsebranch]

test1[x_]:=If[x>10,2x,x^2]

test1[3]

test1[11]

For - klasszikus ciklus (s,i glob.)
For[init,cond,inc,body]

For[s=0;i=0,i≤10,i++,s=s+i];s

Adjuk meg A1,B1 LU felbontását az algoritmussal és ellenőrizzük a felbontás helyességét!

A1={{1,2,3},{4,5,9},{7,8,9}};

MatrixForm[A1]

B1=Table[1,{i,10},{j,10}];

MatrixForm[B1]

(Javaslat)

A1//MatrixForm

{A1L,A1U}=LU0[A1]

MatrixForm[A1L]

MatrixForm[A1U]

A1L.A1U//MatrixForm

B=Table[1,{i,10},{j,10}]

LU0[B]

L, U globális változók

L.U

MatrixForm[L]

MatrixForm[U]

Teszteljük LU0-t olyan 2x2-es mátrixokkal, amelynek nem létezik (egy.) felbontása!
Például, legyen (({{0, 0}, {1, 1}})) (az outputban szereplő x jelzi a hibát) Indoklás: a21=l21*u11+1*0 ⟹ nem létezik ilyen l21

LU0[(({{0, 0}, {1, 1}}))]

Minorok tesztelése Gauss Elimináció, PD mátrixok

Elegendő feltétel az LU felbontás létezésére: Balfelső főminorok nemzérók

Feladat: Ellenőrizzük, hogy az alábbi mátrixok esetén létezik felbontás (és ezekre a mátrixokra a főelem kiv. nélküli parketta alg. meghatározza az egyértelmű felbontást)

V[n_]:=Table[i^j,{i,2,n+1},{j,0,n-1}]

A1=(({{1, 2}, {3, 4}}));A2=(({{1, 2, 3}, {4, 5, 9}, {7, 8, 9}}));A3=V[10];

Javaslat

Det[A2]

Take[A2,2,2]

Minors[A2,2]

Table[Take[A2,j,j],{j,Length[A2]}][[2]]//MatrixForm

Map[Det[#]&,Table[Take[A2,j,j],{j,Length[A2]}]]

AllUpperLeftMinorsNonzeroQ[m_List?MatrixQ]:=And@@Map[(Det[#]≠0)&,Table[Take[m,j,j],{j,Length[m]}]]

AllUpperLeftMinorsNonzeroQ[A3]

AllUpperLeftMinorsNonzeroQ[A1]

AllUpperLeftMinorsNonzeroQ[A2]

PD és PSZD mátrixok

Transzponált,  konjugált transzponált

M1=(({{1, 2}, {3, 4}}));M2=(({{1, 2}, {2, 9}}));M3=(({{1, -2+I, I}, {-3-2I, I, -3}, {12, -3, 3}}));M4=(({{1, -2+I, I}, {-2-I, 19, -3}, {-I, -3, 3}}));

{MatrixForm[Transpose[M1]],MatrixForm[Transpose[M2]]}

{MatrixForm[ConjugateTranspose[M3]],MatrixForm[ConjugateTranspose[M4]]}

HermitianQ[m_List?MatrixQ]:=(m===ConjugateTranspose[m])

HermitianQ[M3]

HermitianQ[M4]

Feladatok PD-Teszt [S/sé]

And@@{True,False,True}

And@@{True,True,True}

PDQ[m_List?MatrixQ]:=And@@Map[(Det[#]>0)&,Table[Take[m,j,j],{j,Length[m]}]]

PDQ[M1]

PDQ[M2]

PSZD-e? (({{1, 1, 1}, {1, 1, 1}, {1, 1, 0}}))

{1,2,-4}.(({{1, 1, 1}, {1, 1, 1}, {1, 1, 0}})).{1,2,-4}

Det[(({{1, 1, 1}, {1, 1, 1}, {1, 1, 0}}))]

(R^H)R felbontás parkettázással (Cholesky) [(LL^T)]

(({{■, ■, ■, ■}, {■, ■, ■, ■}, {■, ■, ■, ■}, {■, ■, □, □}}))

Adjuk meg B Cholesky felbontását.

B=(({{1, 2}, {2, 9}}));

□ White-Box: a parketta algoritmus lépésről lépésre.

B=(R^T)R=(({{(Subscript[r, 11]), 0}, {(Subscript[r, 12]), (Subscript[r, 22])}}))(({{(Subscript[r, 11]), (Subscript[r, 12])}, {0, (Subscript[r, 22])}}))
((Subscript[r, 11])^2)=1 ⟹ (Subscript[r, 11])=1 (Subscript[r, 12])=2 , ((Subscript[r, 22])^2)=5 ⟹(Subscript[r, 22])=(Sqrt[5])

■ Black-Box: A Mathematica built–in függvény a felbontás jobboldali komponensét adja.

R=CholeskyDecomposition[B];

MatrixForm[R]

{MatrixForm[Transpose[R]],MatrixForm[R]}

B=(R^H)R

MatrixForm[Transpose[R].R]

Adjuk meg a Cholesky felbontást.

B=(({{1, 0, -1}, {0, 1, -1}, {-1, -1, 2}}));

Megjegyzés.

Det[B]

CholeskyDecomposition[B]

B nem pozitív defintit, a felbontást nem adja meg a built–in függvény.
Nézzük az összes 2x2-es és 1x1-es főminort.

MapIndexed[Take[#,#2][[1]]&,Minors[B,2]]

MapIndexed[Take[#,#2][[1]]&,Minors[B,1]]

B nem PD, de van neki Cholesky felbontása.

R={{1,0,-1},{0,1,-1},{0,0,0}};MatrixForm[R]

Transpose[R].R//MatrixForm

B//MatrixForm

B=(({{1, -2+I, I}, {-2-I, 19, -3}, {-I, -3, 3}}));

R=CholeskyDecomposition[B];

MatrixForm[R]

ConjugateTranspose[R].R//MatrixForm

MatrixForm[B]

B=(({{1, 1, 1}, {1, 1, 1}, {1, 1, 0}}));

A mátrix szimmetrikus, de nincs neki Cholesky felbontása.
A mátrixnak van negatív főminora!

MapIndexed[Take[#,#2][[1]]&,Minors[B,2]]