and least squares¶
We now look at several decompositions (or factorisations)
of a matrix into products of structured matrices, and their use in solving least squares problems and linear systems.
For a square or rectangular matrix $A ∈ ℝ^{m × n}$ with more rows than columns ($m ≥ n$), we consider:
Copyright By PowCoder代写 加微信 powcoder
The QR decomposition
A = Q R = \underbrace{\begin{bmatrix} 𝐪_1 | \cdots | 𝐪_m \end{bmatrix}}_{m × m} \underbrace{\begin{bmatrix} × & \cdots & × \\ & ⋱ & ⋮ \\ && × \\ &&0 \\ &&⋮ \\ && 0 \end{bmatrix}}_{m × n}
where $Q$ is orthogonal ($Q^⊤Q = I$, $𝐪_j ∈ ℝ^m$) and $R$ is right triangular, which means it
is only nonzero on or to the right of the diagonal.
The reduced QR decomposition
A = Q̂ R̂ = \underbrace{\begin{bmatrix} 𝐪_1 | \cdots | 𝐪_n \end{bmatrix}}_{m × n} \underbrace{\begin{bmatrix} × & \cdots & × \\ & ⋱ & ⋮ \\ && × \end{bmatrix}}_{n × n}
where $Q$ has orthogonal columns ($Q^⊤Q = I$, $𝐪_j ∈ ℝ^m$) and $R̂$ is upper triangular.
For a square matrix we consider the PLU decomposition:
A = P^⊤ LU
where $P$ is a permutation matrix, $L$ is lower triangular and $U$ is upper triangular.
Finally, for a square, symmetric positive definite ($𝐱^⊤ A 𝐱 > 0$ for all $𝐱 ∈ ℝ^n$, $𝐱 \neq 0$)
matrix we consider the Cholesky decomposition:
The importance of these decomposition for square matrices is that their component pieces are easy to invert on a computer:
\begin{align*}
A = P^⊤ LU &\Rightarrow\qquad A^{-1}𝐛 = U^{-1} L^{-1} P 𝐛 \\
A = QR &\Rightarrow\qquad A^{-1}𝐛 = R^{-1} Q^\top 𝐛 \\
A = L L^⊤ &\Rightarrow\qquad A^{-1}𝐛 = L^{-⊤} L^{-1} 𝐛
\end{align*}
and we saw last lecture that triangular and orthogonal matrices are easy to invert when applied to a vector $𝐛$,
e.g., using forward/back-substitution.
For rectangular matrices we will see that they lead to efficient solutions to the least squares problem: find
$𝐱$ that minimizes the 2-norm
\| A 𝐱 – 𝐛 \|.
In this lecture we discuss the followng:
QR and least squares: We discuss the QR decomposition and its usage in solving least squares problems.
Reduced QR and Gram–Schmidt: We discuss computation of the Reduced QR decomposition using Gram–Schmidt.
Householder reflections and QR: We discuss comuting the QR decomposition using Householder reflections.
PLU decomposition: we discuss how the LU decomposition can be computed using Gaussian elimination, and the computation of
the PLU decomposition via Gaussian elimination with pivoting.
Cholesky decomposition: we introduce symmetric positive definite matrices and show that their LU decomposition can be re-interpreted
as a Cholesky decomposition.
Timings: we see the relative trade-off in speed between the different decompositions.
using LinearAlgebra, Plots, BenchmarkTools
1. QR and least squares¶
Here we consider rectangular matrices with more rows than columns.
A QR decomposition decomposes a matrix into an orthogonal matrix $Q$ times a right triangular matrix $R$.
Note the QR decomposition contains within it the reduced QR decomposition:
A = QR = \begin{bmatrix} Q̂ | 𝐪_{n+1} | ⋯ | 𝐪_m \end{bmatrix} \begin{bmatrix} R̂ \\ 𝟎_{m-n × n} \end{bmatrix} = Q̂ R̂.
We can use it to solve a least squares problem using the norm-preserving property (see PS3) of orthogonal matrices:
\| A 𝐱 – 𝐛 \| = \| Q R 𝐱 – 𝐛 \| = \| R 𝐱 – Q^⊤ 𝐛 \| = \left \|
\begin{bmatrix} R̂ \\ 𝟎_{m-n × n} \end{bmatrix} 𝐱 – \begin{bmatrix} Q̂^⊤ \\ 𝐪_{n+1}^⊤ \\ ⋮ \\ 𝐪_m^⊤ \end{bmatrix} 𝐛 \right \|
Now note that the rows $k > n$ are independent of $𝐱$ and are a fixed contribution. Thus to minimise this norm it suffices to
drop them and minimise:
\| R̂ 𝐱 – Q̂^⊤ 𝐛 \|
This norm is minimisable if it is attained. Provided the column rank of $A$ is full, $R̂$ will be invertible (Exercise: why is this?).
Thus we have the solution
𝐱 = R̂^{-1} Q̂^⊤ 𝐛
Example (quadratic fit) Suppose we want to fit noisy data by a quadratic
p(x) = p₀ + p₁ x + p₂ x^2
That is, we want to choose $p₀,p₁,p₂$ at data samples $x_1, \ldots, x_m$ so that the following is true:
p₀ + p₁ x_k + p₂ x_k^2 ≈ f_k
where $f_k$ are given by data. We can reinterpret this as a least squares problem: minimise the norm
\left\| \begin{bmatrix} 1 & x_1 & x_1^2 \\ ⋮ & ⋮ & ⋮ \\ 1 & x_m & x_m^2 \end{bmatrix}
\begin{bmatrix} p₀ \\ p₁ \\ p₂ \end{bmatrix} – \begin{bmatrix} f_1 \\ ⋮ \\ f_m \end{bmatrix} \right \|
We can solve this using the QR decomposition:
m,n = 100,3
x = range(0,1; length=m) # 100 points
f = 2 .+ x .+ 2x.^2 .+ 0.1 .* randn.() # Noisy quadratic
A = x .^ (0:2)’ # 100 x 3 matrix, equivalent to [ones(m) x x.^2]
Q,R̂ = qr(A)
Q̂ = Q[:,1:n] # Q represents full orthogonal matrix so we take first 3 columns
p₀,p₁,p₂ = R̂ \ Q̂’f
3-element Vector{Float64}:
1.991302234012823
1.101309857098941
1.8903486392901003
We can visualise the fit:
p = x -> p₀ + p₁*x + p₂*x^2
scatter(x, f; label=”samples”, legend=:bottomright)
plot!(x, p.(x); label=”quadratic”)
Note that \ with a rectangular system does least squares by default:
3-element Vector{Float64}:
1.9913022340128237
1.1013098570989412
1.8903486392901003
2. Reduced QR and Gram–Schmidt¶
How do we compute the QR decomposition? We begin with a method
you may have seen before in another guise. Write
A = \begin{bmatrix} 𝐚_1 | \dots | 𝐚_n \end{bmatrix}
where $𝐚_k \in ℝ^m$ and assume they are linearly independent ($A$ has full column rank).
Note that the column span of the first $j$ columns of $A$
will be the same as the first $j$ columns of $Q̂$, as
$R̂$ must be non-singular:
\hbox{span}(𝐚_1,\ldots,𝐚_j) = \hbox{span}(𝐪_1,\ldots,𝐪_j)
In other words: the columns of $Q̂$ are an orthogonal basis
of the column span of $A$.
To see this note that since R̂ is triangular we have
\begin{bmatrix} 𝐚_1 | \dots | 𝐚_j \end{bmatrix} = \begin{bmatrix} 𝐪_1 | \dots | 𝐪_j \end{bmatrix} R̂[1:j,1:j]
for all $j$. That is, if $𝐯 ∈ \hbox{span}(𝐚_1,\ldots,𝐚_j)$ we have for $𝐜 ∈ ℝ^j$
𝐯 = \begin{bmatrix} 𝐚_1 | \dots | 𝐚_j \end{bmatrix} 𝐜 = \begin{bmatrix} 𝐪_1 | \dots | 𝐪_j \end{bmatrix} R̂[1:j,1:j] 𝐜 ∈ \hbox{span}(𝐪_1,\ldots,𝐪_j)
while if $𝐰 ∈ \hbox{span}(𝐪_1,\ldots,𝐪_j)$ we have for $𝐝 ∈ ℝ^j$
𝐰 = \begin{bmatrix} 𝐪_1 | \dots | 𝐪_j \end{bmatrix} 𝐝 = \begin{bmatrix} 𝐚_1 | \dots | 𝐚_j \end{bmatrix} R̂[1:j,1:j]^{-1} 𝐝 ∈ \hbox{span}(𝐚_1,\ldots,𝐚_j).
It is possible to find an orthogonal basis using the Gram–Schmidt algorithm,
We construct it via induction:
assume that
\hbox{span}(𝐚_1,\ldots,𝐚_{j-1}) = \hbox{span}(𝐪_1,\ldots,𝐪_{j-1})
where $𝐪_1,\ldots,𝐪_{j-1}$ are orthogonal:
𝐪_k^\top 𝐪_ℓ = δ_{kℓ} = \begin{cases} 1 & k = ℓ \\
0 & \hbox{otherwise} \end{cases}.
for $k,ℓ < j$.
𝐯_j := 𝐚_j - \sum_{k=1}^{j-1} \underbrace{𝐪_k^\top 𝐚_j}_{r_{kj}} 𝐪_k
so that for $k < j$
𝐪_k^\top 𝐯_j = 𝐪_k^\top 𝐚_j - \sum_{k=1}^{j-1} \underbrace{𝐪_k^\top 𝐚_j}_{r_{kj}} 𝐪_k^\top 𝐪_k = 0.
Then we define
𝐪_j := {𝐯_j \over \|𝐯_j\|}.
which sastisfies $𝐪_k^\top 𝐪_j =δ_{kj}$ for $k \leq j$.
We now reinterpret this construction as a reduced QR decomposition.
r_{jj} := {\|𝐯_j\|}
Then rearrange the definition we have
𝐚_j = \begin{bmatrix} 𝐪_1|\cdots|𝐪_j \end{bmatrix} \begin{bmatrix} r_{1j} \\ ⋮ \\ r_{jj} \end{bmatrix}
\begin{bmatrix} 𝐚_1|\cdots|𝐚_j \end{bmatrix}
\begin{bmatrix} r_{11} & \cdots & r_{1j} \\ & ⋱ & ⋮ \\ && r_{jj} \end{bmatrix}
That is, we are computing the reduced QR decomposition column-by-column.
Running this algorithm to $j = n$ completes the decomposition.
Gram–Schmidt in action¶
We are going to compute the reduced QR of a random matrix
A = randn(m,n)
Q,R̂ = qr(A)
Q̂ = Q[:,1:n]
5×4 Matrix{Float64}:
-0.630482 -0.426202 -0.101155 -0.530133
-0.524428 -0.262607 -0.266452 0.482944
-0.129158 0.425953 -0.719869 0.262994
-0.0662478 0.485012 -0.277457 -0.625329
-0.553531 0.576816 0.568838 0.159755
The first column of Q̂ is indeed a normalised first column of A:
R = zeros(n,n)
Q = zeros(m,n)
R[1,1] = norm(A[:,1])
Q[:,1] = A[:,1]/R[1,1]
5-element Vector{Float64}:
0.630482400681896
0.5244283950923772
0.12915824346683272
0.0662478315824829
0.5535306439632841
We now determine the next entries as
R[1,2] = Q[:,1]'A[:,2]
v = A[:,2] - Q[:,1]*R[1,2]
R[2,2] = norm(v)
Q[:,2] = v/R[2,2]
5-element Vector{Float64}:
-0.4262022226633618
-0.2626074204846619
0.4259525739458973
0.48501164628171584
0.5768163623053414
And the third column is then:
R[1,3] = Q[:,1]'A[:,3]
R[2,3] = Q[:,2]'A[:,3]
v = A[:,3] - Q[:,1:2]*R[1:2,3]
R[3,3] = norm(v)
Q[:,3] = v/R[3,3]
5-element Vector{Float64}:
0.1011551707609835
0.26645180209373276
0.7198692579407874
0.27745691988880977
-0.568838270220028
(Note the signs may not necessarily match.)
We can clean this up as a simple algorithm:
function gramschmidt(A)
m,n = size(A)
m ≥ n || error("Not supported")
R = zeros(n,n)
Q = zeros(m,n)
for j = 1:n
for k = 1:j-1
R[k,j] = Q[:,k]'*A[:,j]
v = A[:,j] - Q[:,1:j-1]*R[1:j-1,j]
R[j,j] = norm(v)
Q[:,j] = v/R[j,j]
Q,R = gramschmidt(A)
norm(A - Q*R)
3.521787733787174e-16
Complexity and stability¶
We see within the for j = 1:n loop that we have $O(mj)$ operations. Thus the
total complexity is $O(m n^2)$ operations.
Unfortunately, the Gram–Schmidt algorithm is unstable: the rounding errors when implemented in floating point
accumulate in a way that we lose orthogonality:
A = randn(300,300)
Q,R = gramschmidt(A)
norm(Q'Q-I)
5.455263120264345e-13
3. Householder reflections and QR¶
As an alternative, we will consider using Householder reflections to introduce zeros below
the diagonal.
Thus, if Gram–Schmidt is a process of triangular orthogonalisation (using triangular matrices
to orthogonalise), Householder reflections is a process of orthogonal triangularisation
(using orthogonal matrices to triangularise).
Consider multiplication by the Householder reflection corresponding to the first column,
that is, for
Q_1 := Q_{𝐚_1}^{\rm H},
Q_1 A = \begin{bmatrix} × & × & \cdots & × \\
& × & \cdots & × \\
& ⋮ & ⋱ & ⋮ \\
& × & \cdots & × \end{bmatrix} =
\begin{bmatrix} r_{11} & r_{12} & \cdots & r_{1n} \\
& 𝐚_2^1 & \cdots & 𝐚_n^1 \end{bmatrix}
r_{1j} := (Q_1 𝐚_j)[1] \qquad \hbox{and} \qquad 𝐚_j^1 := (Q_1 𝐚_j)[2:m],
noting $r_{11} = -\hbox{sign}(a_{11})\|𝐚_1\|$ and all entries of $𝐚_1^1$ are zero (thus not included).
That is, we have made the first column triangular.
But now consider
Q_2 := \begin{bmatrix} 1 \\ & Q_{𝐚_2^1}^{\rm H} \end{bmatrix} = Q_{\begin{bmatrix} 0 \\ 𝐚_2^1 \end{bmatrix}}^H
Q_2 Q_1A = \begin{bmatrix} × & × & × & \cdots & × \\
& × & × & \cdots & × \\
&& ⋮ & ⋱ & ⋮ \\
&& × & \cdots & × \end{bmatrix} = \begin{bmatrix} r_{11} & r_{12} & r_{13} & \cdots & r_{1n} \\
& r_{22} & r_{23} & \cdots & r_{2n} \\
&& 𝐚_3^2 & \cdots & 𝐚_n^2 \end{bmatrix}
r_{2j} := (Q_2 𝐚_j^1)[1] \qquad \hbox{and} \qquad 𝐚_j^2 := (Q_2 𝐚_j^1)[2:m-1]
Thus the first two columns are triangular.
The inductive construction is thus clear. If we define $𝐚_j^0 := 𝐚_j$ we
have the construction
\begin{align*}
Q_j &:= \begin{bmatrix} I_{j-1 × j-1} \\ & Q_{𝐚_j^j}^{±,\rm H} \end{bmatrix} \\
𝐚_j^k &:= (Q_k 𝐚_j^{k-1})[2:m-k+1] \\
r_{kj} &:= (Q_k 𝐚_j^{k-1})[1]
\end{align*}
Q_n \cdots Q_1 A = \underbrace{\begin{bmatrix}
r_{11} & \cdots & r_{1n} \\ & ⋱ & ⋮\\
&& r_{nn} \\&& 0 \\ && ⋮ \\ && 0 \end{bmatrix}}_R
A = \underbrace{Q_1 \cdots Q_n}_Q R.
The implementation is cleaner. We do a naive implementation here:
function householderreflection(x)
y = copy(x)
# we cannot use sign(x[1]) in case x[1] == 0
y[1] += (x[1] ≥ 0 ? 1 : -1)*norm(x)
w = y/norm(y)
I - 2*w*w'
function householderqr(A)
m,n = size(A)
R = copy(A)
Q = Matrix(1.0I, m, m)
for j = 1:n
Qⱼ = householderreflection(R[j:end,j])
R[j:end,:] = Qⱼ*R[j:end,:]
Q[:,j:end] = Q[:,j:end]*Qⱼ
A = randn(m, n)
Q,R = householderqr(A)
Note because we are forming a full matrix representation of each Householder
reflection this is a slow algorithm, taking $O(n^4)$ operations. The problem sheet
will consider a better implementation that takes $O(n^3)$ operations.
Example We will now do an example by hand. Consider the $4 × 3$ matrix
A = \begin{bmatrix}
1 & 2 & -1 \\
0 & 15 & 18 \\
-2 & -4 & -4 \\
-2 & -4 & -10
\end{bmatrix}
For the first column we have
Q_1 = I - {1 \over 12} \begin{bmatrix} 4 \\ 0 \\ -2 \\ -2 \end{bmatrix} \begin{bmatrix} 4 & 0 & -2 & -2 \end{bmatrix} =
{1 \over 3} \begin{bmatrix}
-1 & 0 & 2 & 2 \\
0 & 3 & 0 & 0 \\
2 & 0 & 2 & -1 \\
2 & 0 & -1 &2
\end{bmatrix}
Q_1 A = \begin{bmatrix} -3 & -6 & -9 \\
& 15 & 18 \\
& 0 & 0 \\
\end{bmatrix}
In this example the next column is already upper-triangular,
but because of our choice of reflection we will end up swapping the sign, that is
Q_2 = \begin{bmatrix} 1 \\ & -1 \\ && 1 \\ &&& 1 \end{bmatrix}
Q_2 Q_1 A = \begin{bmatrix} -3 & -6 & -9 \\
& -15 & -18 \\
& 0 & 0 \\
\end{bmatrix}
The final reflection is
Q_3 = \begin{bmatrix} I_{2 × 2} \\ & I - \begin{bmatrix} 1 \\ -1 \end{bmatrix} \begin{bmatrix} 1 & -1 \end{bmatrix}
\end{bmatrix} = \begin{bmatrix} å1 \\ & 1 \\ & & 0 & 1 \\ & & 1 & 0 \end{bmatrix}
Q_3 Q_2 Q_1 A = \underbrace{\begin{bmatrix} -3 & -6 & -9 \\
& -15 & -18 \\
& & -6 \\
\end{bmatrix}}_R
A = Q_1 Q_2 Q_3 R = \underbrace{{1 \over 3} \begin{bmatrix}
-1 & 0 & 2 & 2 \\
0 & 3 & 0 & 0 \\
2 & 0 & -1 & 2 \\
2 & 0 & 2 &-1
\end{bmatrix}}_Q \underbrace{\begin{bmatrix} -3 & -6 & -9 \\
& -15 & -18 \\
& & -6 \\
\end{bmatrix}}_R = \underbrace{{1 \over 3} \begin{bmatrix}
-1 & 0 & 2 \\
0 & 3 & 0 \\
2 & 0 & -1 \\
\end{bmatrix}}_Q̂ \underbrace{\begin{bmatrix} -3 & -6 & -9 \\
& -15 & -18 \\
\end{bmatrix}}_R̂
2. PLU Decomposition¶
Just as Gram–Schmidt can be reinterpreted as a reduced QR decomposition,
Gaussian elimination with pivoting can be interpreted as a PLU decomposition.
Special "one-column" lower triangular matrices¶
Consider the following set of $n × n$ lower triangular
matrices which equals identity apart from one-column:
{\cal L}_j := \left\{I + \begin{bmatrix} 𝟎_j \\ 𝐥_j \end{bmatrix} 𝐞_j^⊤ : 𝐥_j ∈ ℝ^{n-j} \right\}
where $𝟎_j$ denotes the zero vector of length $j$.
That is, if $L_j ∈ {\cal L}_j$ then it is equal to the identity matrix apart from in the $j$ th column:
L_j = \begin{bmatrix}
1 \\ & {⋱} \\ && 1 \\
&& ℓ_{j+1,j} & 1 \\
&& ⋮ && \dots \\
&& ℓ_{n,j} & & & 1 \end{bmatrix} =
These satisify the following special properties:
Proposition (one-column lower triangular inverse)
If $L_j \in {\cal L}_j$ then
L_j^{-1} = I - \begin{bmatrix} 𝟎_j \\ 𝐥_j \end{bmatrix} 𝐞_j^⊤ = \begin{bmatrix}
1 \\ & ⋱ \\ && 1 \\
&&-ℓ_{j+1,j} & 1 \\
&&⋮ && \dots \\
&&-ℓ_{n,j} & & & 1 \end{bmatrix} ∈ {\cal L}_j.
Proposition (one-column lower triangular multiplication)
If $L_j \in {\cal L}_j$ and $L_k \in {\cal L}_k$ for $k ≥ j$ then
L_j L_k = I + \begin{bmatrix} 𝟎_j \\ 𝐥_j \end{bmatrix} 𝐞_j^⊤ + \begin{bmatrix} 𝟎_k \\ 𝐥_k \end{bmatrix} 𝐞_k^⊤
Lemma (one-column lower triangular with pivoting)
If $σ$ is a permutation that leaves the first $j$
rows fixed (that is, $σ_ℓ = ℓ$ for $ℓ ≤ j$) and $L_j ∈ {\cal L}_j$ then
P_σ L_j= \tilde L_j P_σ
where $\tilde L_j ∈ {\cal L}_j$.
P_σ = \begin{bmatrix} I_j \\ & P_τ \end{bmatrix}
where $τ$ is the permutation with Cauchy notation
\begin{pmatrix}
1 & \cdots & n-j \\
σ_{j+1}-j & ⋯ & σ_n-j
\end{pmatrix}
Then we have
P_σ L_j = P_σ + \begin{bmatrix} 𝟎_j \\ P_τ 𝐥_j \end{bmatrix} 𝐞_j^⊤ =
\underbrace{(I + \begin{bmatrix} 𝟎_j \\ P_τ 𝐥_j \end{bmatrix} 𝐞_j^⊤)}_{\tilde L_j} P_σ
noting that $𝐞_j^⊤ P_σ = 𝐞_j^⊤$ (as $σ_j = j$).
LU Decomposition¶
Before discussing pivoting, consider standard Gaussian elimation where one row-reduces
to introduce zeros column-by-column. We will mimick the computation of the QR decomposition
to view this as a triangular triangularisation.
Consider the matrix
L_1 = \begin{bmatrix} 1 \\ -{a_{21} \over a_{11}} & 1 \\ ⋮ &&⋱ \\
-{a_{n1} \over a_{11}} &&& 1
\end{bmatrix} = I - \begin{bmatrix} 0 \\ {𝐚_1[2:n] \over 𝐚_1[1]} \end{bmatrix} 𝐞_1^\top.
We then have
L_1 A = \begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\
& 𝐚_2^1 & \cdots & 𝐚_n^1 \end{bmatrix}
where $𝐚_j^1 := (L_1 𝐚_j)[2:n]$ and $u_{1j} = a_{1j}$. But now consider
L_2 := I - \begin{bmatrix} 0 \\ {𝐚_2^1[2:n-1] \over 𝐚_2^1[1]} \end{bmatrix} 𝐞_1^\top.
L_2 L_1 A = \begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\
& u_{22} & u_{23} & \cdots & u_{2n} \\
&& 𝐚_3^2 & \cdots & 𝐚_n^2 \end{bmatrix}
u_{2j} := (𝐚_j^1)[1] \qquad \hbox{and} \qquad 𝐚_j^2 := (L_2 𝐚_j^1)[2:n-1]
Thus the first two columns are triangular.
The inductive construction is again clear. If we define $𝐚_j^0 := 𝐚_j$ we
have the construction
\begin{align*}
L_j &:= I - \begin{bmatrix} 𝟎_j \\ {𝐚_{j+1}^j[2:n-j] \over 𝐚_{j+1}^j[1]} \end{bmatrix} 𝐞_j^⊤ \\
𝐚_j^k &:= (L_k 𝐚_j^{k-1})[2:n-k+1]
u_{kj} &:= (L_k 𝐚_j^{k-1})[1]
\end{align*}
L_{n-1} \cdots L_1 A = \underbrace{\begin{bmatrix}
u_{11} & \cdots & u_{1n} \\ & ⋱ & ⋮\\
&& u_{nn}\end{bmatrix}}_U
A = \underbrace{L_{1}^{-1} \cdots L_{n-1}^{-1}}_L U.
L_j = I + \begin{bmatrix} 𝟎_j \\ \ell_{j+1,j} \\ ⋮ \\ \ell_{n,j} \end{bmatrix} 𝐞_j^\top
and using the properties of inversion and multiplication we therefore deduce
L = \begin{bmatrix} 1 \\ -\ell_{21} & 1 \\
-\ell_{31} & -\ell_{32} & 1 \\
⋮ & ⋮ & ⋱ & ⋱ \\
-\ell_{n1} & -\ell_{n2} & \cdots & -\ell_{n,n-1} & 1
\end{bmatrix}
Example (computer)
We will do a numerical example (by-hand is equivalent to Gaussian elimination).
The first lower triangular matrix is:
A = randn(n,n)
L₁ = I -[0; A[2:end,1]/A[1,1]] * [1 zeros(1,n-1)]
4×4 Matrix{Float64}:
1.0 -0.0 -0.0 -0.0
2.26919 1.0 -0.0 -0.0
-0.42621 -0.0 1.0 -0.0
1.84039 -0.0 -0.0 1.0
Which indeed introduces zeros in the first column:
4×4 Matrix{Float64}:
-0.476424 0.150606 0.472091 0.712055
0.0 -0.410211 1.08783 2.35194
0.0 -1.71671 3.27631 -1.8034
0.0 1.30315 -1.18916 0.937602
Now we make the next lower triangular operator:
L₂ = I - [0; 0; A₁[3:end,2]/A₁[2,2]] * [0 1 zeros(1,n-2)]
4×4 Matrix{Float64}:
1.0 -0.0 -0.0 -0.0
-0.0 1.0 -0.0 -0.0
-0.0 -4.18495 1.0 -0.0
-0.0 3.17677 -0.0 1.0
A₂ = L₂*L₁*A
4×4 Matrix{Float64}:
-0.476424 0.150606 0.472091 0.712055
0.0 -0.410211 1.08783 2.35194
-4.71845e-16 0.0 -1.27622 -11.6462
-8.88178e-16 2.22045e-16 2.26663 8.40917
The last one is:
L₃ = I - [0; 0; 0; A₂[4:end,3]/A₂[3,3]
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com