Exercise: Algorithm to invert matrix
Contents
7. Exercise: Algorithm to invert matrix#
\[
\newcommand{\E}{\text{E}}
\newcommand{\mbf}{\mathbf}
\newcommand{\bs}{\boldsymbol}
\newcommand{\Var}{\text{Var}}
\newcommand{\Cov}{\text{Cov}}
\newcommand{\e}{\frac{1}{\sigma^2_e}}
\newcommand{\f}{\frac{1}{\sigma^2_{\alpha}}}
\]
The formula
\[\begin{split}
\mbf{A}^{-1} =
\begin{bmatrix}
\mbf{A}_{11}^{-1} & \mbf{0} \\
\mbf{0} & \mbf{0}
\end{bmatrix}
+
\begin{bmatrix}
\mbf{q}\mbf{q}' & -\mbf{q} \\
-\mbf{q}' & 1
\end{bmatrix}\frac{1}{r}
\end{split}\]
can also be used to compute \(\mbf{A}^{-1}\) iteratively by starting with \(\mbf{A}_{11} = a_{11}\) and expanding it one row and column at a time.
Consider, for example, the matrix:
\[\begin{split}
\mathbf{A} =
\begin{bmatrix}
2.0 & 0.5 & 0.2 \\
0.5 & 3.0 & 0.1 \\
0.2 & 0.1 & 4.0
\end{bmatrix}
\end{split}\]
A = [2.0 0.5 0.2
0.5 3.0 0.1
0.2 0.1 4.0]
3×3 Matrix{Float64}:
2.0 0.5 0.2
0.5 3.0 0.1
0.2 0.1 4.0
# Inverse of upper 1x1 matrix
A11 = A[1,1]
A11i = 1/A11
0.5
# Inverse of upper 2x2 matrix
A22 = A[2,2]
A12 = A[1,2]
q = A11i*A12
r = A22 - A12'A11i*A12
ri = 1/r
Ai = [A11i .+ ri*q*q' -q*ri
-q'*ri ri]
2×2 Matrix{Float64}:
0.521739 -0.0869565
-0.0869565 0.347826
7.1. Check with Julia inverse function#
s=1:2
inv(A[s,s])
2×2 Matrix{Float64}:
0.521739 -0.0869565
-0.0869565 0.347826
# Inverse of upper 3x3 matrix
s1 = 1:2
s2 = 3
A22 = A[s2,s2]
A12 = A[s1,s2]
A11i = Ai
q = A11i*A12
r = A22 - A12'A11i*A12
ri = 1/r
Ai = [A11i .+ ri*q*q' -q*ri
-q'*ri ri]
3×3 Matrix{Float64}:
0.524038 -0.0865385 -0.0240385
-0.0865385 0.347902 -0.00437063
-0.0240385 -0.00437063 0.251311
7.2. Check with Julia inverse function#
s=1:3
inv(A[s,s])
3×3 Matrix{Float64}:
0.524038 -0.0865385 -0.0240385
-0.0865385 0.347902 -0.00437063
-0.0240385 -0.00437063 0.251311
7.3. Exercise#
Write an iterative function in Julia to invert a symmetric positive definite matrix. Compare the speed of this function with that of the recursive function.
n = 1000
p = 1000
M=rand(n,p)
A=M'M;
function invA(A)
n = size(A,1)
if n==1
return 1/A[1,1]
else
Ai = zeros(n,n) # inverse of A
s1 = 1:n-1
s2 = n
A11 = A[s1,s1]
A11i = invA(A11)
A12 = A[s1,s2]
A22 = A[s2,s2]
q = A11i*A12
r = A22 - q'A11*q
ri = 1/r
Ai = [A11i .+ ri*q*q' -q*ri
-q'*ri ri]
return Ai
end
end
invA (generic function with 1 method)
function invAiter(A)
n = size(A,1)
Ai = 1.0/A[1,1]
for i=2:n
A11i = Ai
s1 = 1:i-1
s2 = i
A11 = A[s1,s1]
A12 = A[s1,s2]
A22 = A[s2,s2]
L = A11i*A12
R = A22 - L'A11*L
Ri = 1/R[1,1]
Ai = [A11i .+ (Ri.* L*L') -Ri.*L
-Ri.*L' Ri.*1]
end
return Ai
end
invAiter (generic function with 1 method)
# complete function
function invAME(A)
n = size(A,1)
Ai = zeros(n,n)
Ai[1,1] = 1.0/A[1,1]
for i=2:n
s1 = 1:i-1
s2 = i
A12 = A[s1,s2]
A22 = A[s2,s2]
L = Ai[s1,s1]*A12
R = A22 - A12'Ai[s1,s1]*A12
Ri = 1/R[1,1]
Ai[s1,s1] += (Ri.* L*L')
Ai[s1,s2] = -Ri.*L
Ai[s2,s1] = -Ri.*L'
Ai[s2,s2] = Ri.*1
end
return Ai
end
invAME (generic function with 1 method)
# complete function
function invAME1(A)
n = size(A,1)
Ai = zeros(n,n)
Ai[1,1] = 1.0/A[1,1]
for i=2:n
s1 = 1:i-1
s2 = i
A12 = A[s1,s2]
A22 = A[s2,s2]
L = Ai[s1,s1]*A12
R = A22 - A12'Ai[s1,s1]*A12
Ri = 1/R[1,1]
for j in s1, k in s1
Ai[j,k] += Ri*L[j]*L[k]
end
Ai[s1,s2] = -Ri.*L
Ai[s2,s1] = -Ri.*L
Ai[s2,s2] = Ri.*1
end
return Ai
end
invAME1 (generic function with 1 method)
@time invA(A);
7.178704 seconds (653.17 k allocations: 14.961 GiB, 12.00% gc time, 1.57% compilation time)
@time invAiter(A);
3.555474 seconds (2.80 M allocations: 12.581 GiB, 16.05% gc time, 10.99% compilation time)
@time invAME(A);
3.374966 seconds (307.23 k allocations: 12.444 GiB, 14.53% gc time, 4.20% compilation time)
@time invAME1(A);
2.111704 seconds (159.73 k allocations: 4.994 GiB, 10.39% gc time, 3.97% compilation time)
@time inv(A);
0.068327 seconds (5 allocations: 8.126 MiB)
a = [1.0 2.0
3.0 4.0]
2×2 Matrix{Float64}:
1.0 2.0
3.0 4.0
v = [5.0,6.0]
2-element Vector{Float64}:
5.0
6.0
a[2,:] = v
2-element Vector{Float64}:
5.0
6.0
a
2×2 Matrix{Float64}:
1.0 2.0
5.0 6.0