CS代考 and least squares¶

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