462 lines
13 KiB
Plaintext
462 lines
13 KiB
Plaintext
---
|
||
engine: julia
|
||
---
|
||
|
||
```{julia}
|
||
#| error: false
|
||
#| echo: false
|
||
#| output: false
|
||
using InteractiveUtils
|
||
```
|
||
|
||
# Linear Algebra in Julia
|
||
|
||
```{julia}
|
||
using LinearAlgebra
|
||
```
|
||
|
||
The `LinearAlgebra` package provides among other things:
|
||
|
||
- additional subtypes of `AbstractMatrix`: usable like other matrices, e.g.,
|
||
|
||
- `Tridiagonal`
|
||
- `SymTridiagonal`
|
||
- `Symmetric`
|
||
- `UpperTriangular`
|
||
|
||
|
||
- additional/extended functions: `norm`, `opnorm`, `cond`, `inv`, `det`, `exp`, `tr`, `dot`, `cross`, ...
|
||
|
||
- a universal solver for systems of linear equations: `\`
|
||
- `x = A \ b` solves $A \mathbf{x}=\mathbf{b}$ by appropriate matrix factorization and forward/backward substitution
|
||
|
||
- [Matrix factorizations](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#man-linalg-factorizations)
|
||
- `LU`
|
||
- `QR`
|
||
- `Cholesky`
|
||
- `SVD`
|
||
- ...
|
||
|
||
- Computation of eigenvalues/eigenvectors
|
||
|
||
- `eigen`, `eigvals`, `eigvecs`
|
||
|
||
|
||
- Access to BLAS/LAPACK functions
|
||
|
||
## Matrix Types
|
||
|
||
|
||
|
||
```{julia}
|
||
A = SymTridiagonal(fill(1.0, 4), fill(-0.3, 3))
|
||
```
|
||
|
||
|
||
```{julia}
|
||
B = UpperTriangular(A)
|
||
```
|
||
|
||
These types are stored space-efficiently. The usual arithmetic operations are implemented:
|
||
|
||
```{julia}
|
||
A + B
|
||
```
|
||
|
||
Read-only index access is possible,
|
||
```{julia}
|
||
A[1,4]
|
||
```
|
||
|
||
write operations not necessarily:
|
||
```{julia}
|
||
#| error: true
|
||
A[1,3] = 17
|
||
```
|
||
|
||
Conversion to a 'normal' matrix is possible, for example, with `collect()`:
|
||
```{julia}
|
||
A2 = collect(A)
|
||
```
|
||
|
||
### The Identity Matrix `I`
|
||
|
||
`I` denotes an identity matrix (square, diagonal elements = 1, all others = 0) of the respective required size
|
||
|
||
|
||
```{julia}
|
||
A + 4I
|
||
```
|
||
|
||
|
||
## Norms
|
||
|
||
To be able to study questions such as conditioning or convergence of an algorithm, we need a metric. For linear spaces, it is appropriate to define the metric via a norm:
|
||
$$
|
||
d(x,y) := ||x-y||
|
||
$$
|
||
|
||
### $p$-Norms
|
||
|
||
A simple class of norms in $ℝ^n$ are the $p$-norms
|
||
$$
|
||
||\mathbf{x}||_p = \left(\sum |x_i|^p\right)^\frac{1}{p},
|
||
$$
|
||
which generalize the Euclidean norm $p=2$.
|
||
|
||
|
||
:::{.callout-note}
|
||
|
||
## The Max-Norm $p=\infty$
|
||
|
||
Let $x_{\text{max}}$ be the _largest in absolute value_ component of $\mathbf{x}\in ℝ^n$. Then always
|
||
$$ |x_{\text{max}}| \le ||\mathbf{x}||_p \le n^\frac{1}{p} |x_{\text{max}}|
|
||
$$
|
||
(Consider a vector whose components are all equal to $x_{\text{max}}$ respectively a vector whose components are all equal to zero except $x_{\text{max}}$.)
|
||
|
||
Thus follows
|
||
$$
|
||
\lim_{p \rightarrow \infty} ||\mathbf{x}||_p = |x_{\text{max}}| =: ||\mathbf{x}||_\infty.
|
||
$$
|
||
:::
|
||
|
||
|
||
In Julia, the `LinearAlgebra` package defines a function `norm(v, p)`.
|
||
|
||
```{julia}
|
||
v = [3, 4]
|
||
w = [-1, 2, 33.2]
|
||
|
||
@show norm(v) norm(v, 2) norm(v, 1) norm(v, 4) norm(w, Inf);
|
||
```
|
||
|
||
- If the 2nd argument `p` is missing, `p=2` is set.
|
||
- The 2nd argument can also be `Inf` (i.e., $+\infty$).
|
||
- The 1st argument can be any container full of numbers. The sum $\sum |x_i|^p$ extends over *all* elements of the container.
|
||
- Thus, for a matrix `norm(A)` is equal to the _Frobenius norm_ of the matrix `A`.
|
||
|
||
```{julia}
|
||
A = [1 2 3
|
||
4 5 6
|
||
7 8 9]
|
||
norm(A) # Frobenius norm
|
||
```
|
||
|
||
|
||
|
||
Since norms are homogeneous under multiplication with scalars,
|
||
$||\lambda \mathbf{x}|| = |\lambda|\cdot||\mathbf{x}||$, they are completely determined by the specification of the unit ball. Subadditivity of the norm (triangle inequality) is equivalent to the convexity of the unit ball
|
||
(code visible by clicking).
|
||
```{julia}
|
||
#| code-fold: true
|
||
#| fig-cap: "Unit balls in $ℝ^2$ for different $p$-norms: $p$=0.8; 1; 1.5; 2; 3.001 and 1000"
|
||
using Plots
|
||
|
||
colors=[:purple, :green, :red, :blue,:aqua, :black]
|
||
x=LinRange(-1, 1, 1000)
|
||
y=LinRange(-1, 1, 1000)
|
||
|
||
fig1=plot()
|
||
for p ∈ (0.8, 1, 1.5, 2, 3.001, 1000)
|
||
contour!(x,y, (x,y) -> p * norm([x, y], p), levels=[p], aspect_ratio=1,
|
||
cbar=false, color=[pop!(colors)], contour_labels=true, ylim=[-1.1, 1.1])
|
||
end
|
||
fig1
|
||
```
|
||
As one can see, $p\ge 1$ must hold so that the unit ball is convex and $||.||_p$ is a norm.
|
||
|
||
However, the Julia function `norm(v, p)` returns a result for arbitrary parameters `p`.
|
||
|
||
|
||
### Induced Norms (Operator Norms)
|
||
|
||
Matrices $A$ represent linear mappings $\mathbf{v}\mapsto A\mathbf{v}$. The matrix norm induced by a vector norm answers the question:
|
||
|
||
> _"By what factor can a vector be maximally stretched by the transformation $A$?"_
|
||
|
||
Due to the homogeneity of the norm under multiplication with scalars, it is sufficient to consider the image of the unit ball under the transformation $A$.
|
||
|
||
::: {.callout-tip}
|
||
## Definition
|
||
Let $V$ be a vector space with dimension $0<n<\infty$ and
|
||
$A$ an $n\times n$-matrix. Then
|
||
$$
|
||
||A||_p = \max_{||\mathbf{v}||_p=1} ||A\mathbf{v}||_p
|
||
$$
|
||
:::
|
||
|
||
Induced norms can only be calculated with difficulty for general $p$. Exceptions are the cases
|
||
|
||
- $p=1$: column sum norm
|
||
- $p=2$: spectral norm and
|
||
- $p=\infty$: row sum norm
|
||
|
||
These 3 cases are implemented in Julia in the function `opnorm(A, p)` from the `LinearAlgebra` package, where again `opnorm(A) = opnorm(A, 2)` holds.
|
||
|
||
```{julia}
|
||
A = [ 0 1
|
||
1.2 1.5 ]
|
||
|
||
@show opnorm(A, 1) opnorm(A, Inf) opnorm(A, 2) opnorm(A);
|
||
```
|
||
|
||
The following picture shows the effect of $A$ on unit vectors. Vectors of the same color are mapped onto each other. (Code visible by clicking):
|
||
|
||
```{julia}
|
||
#| error: false
|
||
#| echo: false
|
||
#| output: false
|
||
|
||
using CairoMakie
|
||
CairoMakie.activate!(type = "svg")
|
||
# set_theme!(size=(1600, 360))
|
||
```
|
||
|
||
```{julia}
|
||
#| code-fold: true
|
||
#| fig-cap: "Image of the unit ball under $v \\mapsto Av$ with $||A||\\approx 2.088$"
|
||
|
||
using CairoMakie
|
||
|
||
# Makie bug https://github.com/MakieOrg/Makie.jl/issues/3255
|
||
# Workaround https://github.com/MakieOrg/Makie.jl/issues/2607#issuecomment-1385816645
|
||
tri = BezierPath([
|
||
MoveTo(Point2f(-0.5, -1)), LineTo(0, 0), LineTo(0.5, -1), ClosePath()
|
||
])
|
||
|
||
A = [ 0 1
|
||
1.2 1.5 ]
|
||
|
||
t = LinRange(0, 1, 100)
|
||
xs = sin.(2π*t)
|
||
ys = cos.(2π*t)
|
||
Axys = A * [xs, ys]
|
||
|
||
u = [sin(n*π/6) for n=0:11]
|
||
v = [cos(n*π/6) for n=0:11]
|
||
x = y = zeros(12)
|
||
|
||
Auv = A * [u,v]
|
||
|
||
fig2 = Figure(size=(800, 400))
|
||
lines(fig2[1, 1], xs, ys, color=t, linewidth=5, colormap=:hsv, axis=(; aspect = 1, limits=(-2,2, -2,2),
|
||
title=L"$\mathbf{v}$", titlesize=30))
|
||
arrows!(fig2[1,1], x, y, u, v, arrowsize=10, arrowhead=tri, colormap=:hsv, linecolor=range(0,11), linewidth=3)
|
||
|
||
Legend(fig2[1,2], MarkerElement[], String[], L"⟹", width=40, height=30, titlesize=30, framevisible=false)
|
||
|
||
lines(fig2[1,3], Axys[1], Axys[2], color=t, linewidth=5, colormap=:hsv, axis=(; aspect=1, limits=(-2,2, -2,2),
|
||
title=L"$A\mathbf{v}$", titlesize=30))
|
||
arrows!(fig2[1,3], x, y, Auv[1], Auv[2], arrowsize=10, arrowhead=tri, colormap=:hsv, linecolor=range(0,11),
|
||
linewidth=3)
|
||
|
||
fig2
|
||
```
|
||
|
||
### Condition Number
|
||
|
||
For p = 1, p = 2 (default) or p = Inf, `cond(A,p)` returns the condition number in the $p$-norm
|
||
$$
|
||
\text{cond}_p(A) = ||A||_p \cdot ||A^{-1}||_p
|
||
$$
|
||
|
||
```{julia}
|
||
@show cond(A, 1) cond(A, 2) cond(A) cond(A, Inf);
|
||
```
|
||
|
||
## Matrix Factorizations
|
||
|
||
Basic tasks of numerical linear algebra:
|
||
|
||
- Solve a system of linear equations $A\mathbf{x} = \mathbf{b}$.
|
||
- If no solution exists, find the best approximation, i.e., the vector $\mathbf{x}$ that minimizes $||A\mathbf{x} - \mathbf{b}||$.
|
||
- Find eigenvalues and eigenvectors $A\mathbf{x} = \lambda \mathbf{x}$ of $A$.
|
||
|
||
These tasks can be solved using matrix factorizations. Some fundamental matrix factorizations:
|
||
|
||
- **LU decomposition** $A=L\cdot U$
|
||
- factorizes a matrix as a product of a _lower_ and an _upper_ triangular matrix
|
||
- in German also called LR decomposition (but the Julia function is called `lu()`)
|
||
- always works (possibly after row exchanges - pivoting)
|
||
- **Cholesky decomposition** $A=L\cdot L^*$
|
||
- the upper triangular matrix is the conjugate of the lower,
|
||
- half the effort compared to LU
|
||
- only works if $A$ is Hermitian and positive definite
|
||
- **QR decomposition** $A=Q\cdot R$
|
||
- decomposes $A$ as a product of an orthogonal (or unitary in the complex case) matrix and an upper triangular matrix
|
||
- $Q$ is length-preserving (rotations and/or reflections); the scalings are described by $R$
|
||
- always works
|
||
- **SVD** _(Singular value decomposition)_: $A = U\cdot D \cdot V^*$
|
||
- $U$ and $V$ are orthogonal (or unitary), $D$ is a diagonal matrix with entries on the diagonal $σ_i\ge 0$, the so-called _singular values_ of $A$.
|
||
- Every linear transformation $\mathbf{v} \mapsto A\mathbf{v}$ can thus be represented as a rotation (and/or reflection) $V^*$, followed by a pure scaling $v_i \mapsto \sigma_i v_i$ and another rotation $U$.
|
||
|
||
|
||
### LU Factorization
|
||
|
||
LU factorization is Gaussian elimination. The result of Gaussian elimination is the upper triangular matrix $U$. The lower triangular matrix $L$ contains ones on the diagonal and the non-diagonal entries $l_{ij}$ are equal to minus the coefficients by which row $Z_j$ is multiplied and added to row $Z_i$ in the Gaussian algorithm.
|
||
An example:
|
||
$$
|
||
A=\left[
|
||
\begin{array}{ccc}
|
||
1 &2 &2 \\
|
||
3 &-4& 4 \\
|
||
-2 & 1 & 5
|
||
\end{array}\right]
|
||
~ \begin{array}{c}
|
||
~\\
|
||
Z_2 \mapsto Z_2 \mathbin{\color{red}-}\textcolor{red}{3} Z_1\\
|
||
Z_3 \mapsto Z_3 + \textcolor{red}{2} Z_1
|
||
\end{array} \quad \Longrightarrow\
|
||
\left[
|
||
\begin{array}{ccc}
|
||
1 &2 &2 \\
|
||
&-10& -2 \\
|
||
& 5 & 9
|
||
\end{array}\right]
|
||
~ \begin{array}{c}
|
||
~\\
|
||
~\\
|
||
Z_3 \mapsto Z_3 + \textcolor{red}{\frac{1}{2}} Z_2
|
||
\end{array} \quad \Longrightarrow\
|
||
\left[
|
||
\begin{array}{ccc}
|
||
1 &2 &2 \\
|
||
&-10& -2 \\
|
||
& & 8
|
||
\end{array}\right] = U
|
||
$$
|
||
$$
|
||
A = \left[
|
||
\begin{array}{ccc}
|
||
1 && \\
|
||
\mathbin{\color{red}+}\textcolor{red}{3} &1 & \\
|
||
\mathbin{\color{red}-}\textcolor{red}{2} & \mathbin{\color{red}-}\textcolor{red}{\frac{1}{2}}& 1
|
||
\end{array}
|
||
\right] \cdot
|
||
\left[
|
||
\begin{array}{ccc}
|
||
1 &2 &2 \\
|
||
&-10& -2 \\
|
||
& & 8
|
||
\end{array}\right]
|
||
$$
|
||
|
||
|
||
- Often in practice: $A\mathbf{x}=\mathbf{b}$ must be solved for one $A$ and many right-hand sides $\mathbf{b}$.
|
||
- The factorization, whose effort grows cubically $\sim n^3$ with the matrix size $n$, only needs to be done once.
|
||
- The subsequent effort of forward/backward substitution for each $\mathbf{b}$ is only quadratically $\sim n^2$.
|
||
|
||
The `LinearAlgebra` package of Julia contains the function `lu(A, options)` for calculating an LU decomposition:
|
||
```{julia}
|
||
A = [ 1 2 2
|
||
3 -4 4
|
||
-2 1 5]
|
||
|
||
L, U, _ = lu(A, NoPivot())
|
||
display(L)
|
||
display(U)
|
||
```
|
||
|
||
|
||
|
||
#### Pivoting
|
||
|
||
Let's look at one step of Gaussian elimination:
|
||
$$
|
||
\left[
|
||
\begin{array}{cccccc}
|
||
* &\cdots & * & * & \cdots & * \\
|
||
&\ddots & \vdots &\vdots && \vdots \\
|
||
&& * & * &\cdots & * \\
|
||
&& & \textcolor{red}{a_{ij}}&\cdots & a_{in}\\
|
||
&& & \textcolor{blue}{a_{i+1,j}}&\cdots & a_{i+1,n}\\
|
||
&&& \textcolor{blue}{\vdots} &&\vdots \\
|
||
&& & \textcolor{blue}{a_{mj}}&\cdots & a_{mn}
|
||
\end{array}
|
||
\right]
|
||
$$
|
||
The goal is to make the entry $a_{i+1,j}$ disappear by adding an appropriate multiple of row $Z_i$ to row $Z_{i+1}$. This only works if the _pivot element_ $\textcolor{red}{a_{ij}}$ is not zero. If $\textcolor{red}{a_{ij}}=0$, we must exchange rows to fix this.
|
||
|
||
Furthermore, the conditioning of the algorithm is best if we arrange the matrix at each step so that the pivot element is the largest in absolute value in the corresponding column of the remaining submatrix. In (row) pivoting, rows are exchanged at each step to ensure that
|
||
$$
|
||
|\textcolor{red}{a_{ij}}|=\max_{k=i,...,m} |\textcolor{blue}{a_{kj}}|
|
||
$$
|
||
|
||
#### LU in Julia
|
||
|
||
- The factorizations in Julia return a special object that contains the matrix factors and additional information.
|
||
- The Julia function `lu(A)` performs an LU factorization with pivoting.
|
||
|
||
|
||
```{julia}
|
||
F = lu(A)
|
||
typeof(F)
|
||
```
|
||
Elements of the object:
|
||
```{julia}
|
||
@show F.L F.U F.p;
|
||
```
|
||
One can also use an appropriate tuple on the left side:
|
||
```{julia}
|
||
L, U, p = lu(A);
|
||
p
|
||
```
|
||
|
||
The permutation vector indicates how the rows of the matrix have been permuted. It holds: $$ L\cdot U = PA$$. The syntax of indirect indexing allows applying the row permutation with the notation `A[p,:]`:
|
||
```{julia}
|
||
display(A)
|
||
display(A[p,:])
|
||
display(L*U)
|
||
```
|
||
|
||
|
||
Forward/backward substitution with an `LU`-object is performed by the `\` operator:
|
||
```{julia}
|
||
b = [1, 2, 3]
|
||
x = F \ b
|
||
```
|
||
Verification:
|
||
```{julia}
|
||
A * x - b
|
||
```
|
||
|
||
In Julia, the `\` operator hides a quite universal "matrix solver", and it can also be applied directly:
|
||
```{julia}
|
||
A \ b
|
||
```
|
||
An appropriate factorization is performed implicitly, but its result is not saved.
|
||
|
||
|
||
### QR Decomposition
|
||
|
||
|
||
The function `qr()` returns a special QR object that contains the components $Q$ and $R$. The orthogonal (or unitary) matrix $Q$ is
|
||
[in an optimized form](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#man-linalg-abstractq) stored. Conversion to a "normal" matrix is, as always, possible with `collect()` if needed.
|
||
```{julia}
|
||
F = qr(A)
|
||
@show typeof(F) typeof(F.Q)
|
||
display(collect(F.Q))
|
||
display(F.R)
|
||
```
|
||
|
||
### Appropriate Factorization
|
||
|
||
|
||
The function `factorize()` returns a factorization form adapted to the matrix type, see [documentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.factorize) for details.
|
||
If solutions for several right-hand sides $\mathbf{b_1}, \mathbf{b_2},...$ are needed, the factorization should only be performed once:
|
||
|
||
|
||
```{julia}
|
||
Af = factorize(A)
|
||
```
|
||
|
||
|
||
```{julia}
|
||
Af \ [1, 2, 3]
|
||
```
|
||
|
||
|
||
```{julia}
|
||
Af \ [5, 7, 9]
|
||
```
|