## Gauss Jordan, LU¶

### A Bit More about Blocks -- Consider 1x1 by 1x2¶

In [195]:
2*[3 4]

Out[195]:
1x2 Array{Int64,2}:
6  8

In [196]:
[2*3 2*4]

Out[196]:
1x2 Array{Int64,2}:
6  8


### Blocks work the same way¶

In [197]:
A=rand(1:3,2,4); B=rand(1:3,4,3); C=rand(1:3,4,2);
println( A*[B C])
println( [A*B A*C])

15	16	13	15	16
15	16	13	17	14

15	16	13	15	16
15	16	13	17	14


Question: What dimensions of A and B and C allow this to work There are a number of ways to check that with the right dimensions A[B C]=[AB AC] 1. The highest level way is the block view Block view = if blocks conform then block matmul always works like elementwise matmul 2. Also worthwhile to consider the column view The columns of A[B C] are all the matrix-vector products of A*the columns of [B C] Why is this the same as [AB AC], the columns of AB concatanated with the columns of AC? 3. You might try other views of matmul too You might think some are easier on the brain than others, and also some views while less easy on the brain, are more powerful. You get to decide what you like best.

### Gauss Jordan in high level block form M[A I]=[I M] so M is a left inverse (algorithm hidden!)¶

In [198]:
rand3(n)=rand(1:3,n,n);  # define rand3 (A little shorthand)

Out[198]:
# methods for generic function rand3
rand3(n) at In[198]:1

In [200]:
I=eye(3); A=rand3(3); M=inv(A);     # Might be singular, careful (single exception or large numbers indicate singularity)

Out[200]:
3x3 Array{Float64,2}:
0.4  -0.2   0.0
-0.6  -0.2   1.0
0.4   0.8  -1.0

In [209]:
I=eye(3);
showcompact(
M*[A I]    );  println('\n')
showcompact(
[ I M]
)

3x6 Array{Float64,2}:
1.0  0.0  0.0   0.4  -0.2   0.0
0.0  1.0  0.0  -0.6  -0.2   1.0
0.0  0.0  1.0   0.4   0.8  -1.0

3x6 Array{Float64,2}:
1.0  0.0  0.0   0.4  -0.2   0.0
0.0  1.0  0.0  -0.6  -0.2   1.0
0.0  0.0  1.0   0.4   0.8  -1.0


### Gauss Jordan details -- do this with pencil and paper once or twice¶

In [210]:
Z=[ A I]

Out[210]:
3x6 Array{Float64,2}:
3.0  1.0  1.0  1.0  0.0  0.0
1.0  2.0  2.0  0.0  1.0  0.0
2.0  2.0  1.0  0.0  0.0  1.0

In [211]:
Z[2,:] -= (Z[2,1]/Z[1,1])*Z[1,:];Z    # Method 1. (Easier for Programming) Elimination via row operations

Out[211]:
3x6 Array{Float64,2}:
3.0  1.0      1.0       1.0       0.0  0.0
0.0  1.66667  1.66667  -0.333333  1.0  0.0
2.0  2.0      1.0       0.0       0.0  1.0

In [211]:
Z=[ A I];

In [212]:
E=copy(I); E[2,1]=-Z[2,1]/Z[1,1];E  # Method 2: Create an elimination matrix with the negative multiplier in the (2,1) position, say

Out[212]:
3x3 Array{Float64,2}:
1.0       0.0  0.0
-0.333333  1.0  0.0
0.0       0.0  1.0

In [213]:
Z=E*Z

Out[213]:
3x6 Array{Float64,2}:
3.0  1.0      1.0       1.0       0.0  0.0
0.0  1.66667  1.66667  -0.333333  1.0  0.0
2.0  2.0      1.0       0.0       0.0  1.0

##### Continue with method 1¶
In [214]:
Z[3,:] -= (Z[3,1]/Z[1,1])*Z[1,:];Z

Out[214]:
3x6 Array{Float64,2}:
3.0  1.0      1.0        1.0       0.0  0.0
0.0  1.66667  1.66667   -0.333333  1.0  0.0
0.0  1.33333  0.333333  -0.666667  0.0  1.0

##### Column 2¶
In [215]:
Z[1,:] -= (Z[1,2]/Z[2,2])*Z[2,:];Z

Out[215]:
3x6 Array{Float64,2}:
3.0  0.0      0.0        1.2       -0.6  0.0
0.0  1.66667  1.66667   -0.333333   1.0  0.0
0.0  1.33333  0.333333  -0.666667   0.0  1.0

In [216]:
Z[3,:] -= (Z[3,2]/Z[2,2])*Z[2,:];Z

Out[216]:
3x6 Array{Float64,2}:
3.0  0.0       0.0       1.2       -0.6  0.0
0.0  1.66667   1.66667  -0.333333   1.0  0.0
0.0  0.0      -1.0      -0.4       -0.8  1.0

##### Column 3¶
In [217]:
Z[1,:] -= (Z[1,3]/Z[3,3])*Z[3,:];Z

Out[217]:
3x6 Array{Float64,2}:
3.0  0.0       0.0       1.2       -0.6  0.0
0.0  1.66667   1.66667  -0.333333   1.0  0.0
0.0  0.0      -1.0      -0.4       -0.8  1.0

In [218]:
Z[2,:] -= (Z[2,3]/Z[3,3])*Z[3,:];Z

Out[218]:
3x6 Array{Float64,2}:
3.0  0.0       0.0   1.2  -0.6       0.0
0.0  1.66667   0.0  -1.0  -0.333333  1.66667
0.0  0.0      -1.0  -0.4  -0.8       1.0

##### Divide out the pivots¶
In [219]:
Z[1,:]/=Z[1,1];Z

Out[219]:
3x6 Array{Float64,2}:
1.0  0.0       0.0   0.4  -0.2       0.0
0.0  1.66667   0.0  -1.0  -0.333333  1.66667
0.0  0.0      -1.0  -0.4  -0.8       1.0

In [220]:
Z[2,:]/=Z[2,2];Z

Out[220]:
3x6 Array{Float64,2}:
1.0  0.0   0.0   0.4  -0.2  0.0
0.0  1.0   0.0  -0.6  -0.2  1.0
0.0  0.0  -1.0  -0.4  -0.8  1.0

In [221]:
Z[3,:]/=Z[3,3];Z

Out[221]:
3x6 Array{Float64,2}:
1.0   0.0  0.0   0.4  -0.2   0.0
0.0   1.0  0.0  -0.6  -0.2   1.0
-0.0  -0.0  1.0   0.4   0.8  -1.0

In [222]:
inv(A)

Out[222]:
3x3 Array{Float64,2}:
0.4  -0.2   0.0
-0.6  -0.2   1.0
0.4   0.8  -1.0

##### We dit it!¶
In [223]:
Point of view: we solved Ax =[e1 e2 e3] simultaneously giving a right inverse but the block view shows us it is also a left inverse hence a left invese is a right inverse

## Elimination = Factorization = LU¶

In [276]:
A=[ 6.0 5 10; 8 13 8; 3 14 7]; Aoriginal=copy(A);

In [262]:
A[2,:]-= (4/3)*A[1,:]
A

Out[262]:
3x3 Array{Float64,2}:
6.0   5.0      10.0
0.0   6.33333  -5.33333
3.0  14.0       7.0

In [263]:
A[3,:] -= ( 0.5 )*A[1,:]
A

Out[263]:
3x3 Array{Float64,2}:
6.0   5.0      10.0
0.0   6.33333  -5.33333
0.0  11.5       2.0

In [264]:
A[3,:]-=  (11.5/6.333333   )*A[2,:];A

Out[264]:
3x3 Array{Float64,2}:
6.0   5.0         10.0
0.0   6.33333     -5.33333
0.0  -6.05263e-7  11.6842

In [265]:
U=ans

Out[265]:
3x3 Array{Float64,2}:
6.0   5.0         10.0
0.0   6.33333     -5.33333
0.0  -6.05263e-7  11.6842

In [266]:
(11.5/6.333333 )

Out[266]:
1.8157895692520827

In [267]:
69/38

Out[267]:
1.8157894736842106

In [268]:
U

Out[268]:
3x3 Array{Float64,2}:
6.0   5.0         10.0
0.0   6.33333     -5.33333
0.0  -6.05263e-7  11.6842

In [269]:
L=eye(3)

Out[269]:
3x3 Array{Float64,2}:
1.0  0.0  0.0
0.0  1.0  0.0
0.0  0.0  1.0

In [270]:
L[2,1]=4/3; L[3,1]=0.5; L[3,2]=69/38;L

Out[270]:
3x3 Array{Float64,2}:
1.0      0.0      0.0
1.33333  1.0      0.0
0.5      1.81579  1.0

In [271]:
L*U

Out[271]:
3x3 Array{Float64,2}:
6.0   5.0  10.0
8.0  13.0   8.0
3.0  14.0   7.0

In [277]:
Aoriginal

Out[277]:
3x3 Array{Float64,2}:
6.0   5.0  10.0
8.0  13.0   8.0
3.0  14.0   7.0

In [258]:
E

Out[258]:
3x3 Array{Float64,2}:
1.0       0.0  0.0
-0.333333  1.0  0.0
0.0       0.0  1.0

In [237]:
inv(E)

Out[237]:
3x3 Array{Float64,2}:
1.0        0.0  0.0
0.333333   1.0  0.0
-0.0       -0.0  1.0

In [238]:
M=rand(3,3)

Out[238]:
3x3 Array{Float64,2}:
0.0385084  0.992836   0.113299
0.399386   0.757624   0.419412
0.428451   0.0878158  0.299829

In [239]:
E*M

Out[239]:
3x3 Array{Float64,2}:
0.0385084  0.992836   0.113299
0.38655    0.426679   0.381646
0.428451   0.0878158  0.299829

In [240]:
inv(E)* (E*M)

Out[240]:
3x3 Array{Float64,2}:
0.0385084  0.992836   0.113299
0.399386   0.757624   0.419412
0.428451   0.0878158  0.299829

In []: