xxxxxxxxxx
 
begin 
    using Pkg
    Pkg.activate(mktempdir()) 
    Pkg.add("PyPlot") 
    Pkg.add("BenchmarkTools")
    Pkg.add("PlutoUI")
    using PlutoUI,PyPlot,BenchmarkTools
end

6.4 s

xxxxxxxxxx
 
using LinearAlgebra

189 μs

Linear System Solution

Let $A$ be an $N \times N$ matrix, $b \in R^{N}$ .
Solve $A x = b$

Direct methods:

Exact up to machine precision
Sometimes expensive, sometimes not

Iterative methods:

"Only" approximate
With good convergence and proper accuracy control, results may be not worse than for direct methods
Sometimes expensive, sometimes not
Convergence guarantee is problem dependent and can be tricky

13.3 μs

Matrix & vector norms

Vector norms

let $x = (x_{i}) \in R^{n}$

4.5 μs

2.7 μs

$| | x | |_{1} = \sum_{i = 1}^{n} | x_{i} |$ : sum norm, $l_{1}$ -norm

4.2 μs

1.5

xxxxxxxxxx
 
norm(v,1)

75.3 ms

$| | x | |_{2} = \sqrt{\sum_{i = 1}^{n} x_{i}^{2}}$ : Euclidean norm, $l_{2}$ -norm

3.9 μs

1.06301

1.06301

xxxxxxxxxx
 
norm(v,2),norm(v)

2.2 ms

$| | x | |_{\infty} = {max}_{i = 1}^{n} | x_{i} |$ : maximum norm, $l_{\infty}$ -norm

4.0 μs

1.0

xxxxxxxxxx
 
norm(v,Inf)

91.0 ms

Matrix norms

Matrix $A = (a_{i j}) \in R^{n} \times R^{n}$

Representation of linear operator $A : R^{n} \to R^{n}$ defined by $A : x \mapsto y = A x$ with $y_{i} = \sum_{j = 1}^{n} a_{i j} x_{j}$
Vector norm $| | \cdot | |_{p}$ induces corresponding matrix norm:

$| | A | |_{p} = max_{x \in R^{n}, x \neq 0} \frac{| | A x | |_{p}}{| | x | |_{p}} = max_{x \in R^{n}, | | x | |_{p} = 1} \frac{| | A x | |_{p}}{| | x | |_{p}}$

6.5 μs

3×3 Array{Float64,2}:
 3.0  2.0  3.0
 0.1  0.3  0.5
 0.6  2.0  3.0

xxxxxxxxxx
 
M=[3.0 2.0 3.0; 
   0.1 0.3 0.5; 
   0.6 2.0 3.0]

16.4 ms

$| | A | |_{1} = {max}_{j = 1}^{n} \sum_{i = 1}^{n} | a_{i j} |$ maximum of column sums of absolute values of entries

3.2 μs

6.5

xxxxxxxxxx
 
opnorm(M,1)

94.0 ms

$| | A | |_{2} = \sqrt{λ_{m a x}}$ with $λ_{m a x}$ : largest eigenvalue of $A^{T} A$ .

xxxxxxxxxx
 
md"""
- ``||A||_2=\sqrt{\lambda_{max}}`` with ``\lambda_{max}``: largest eigenvalue of ``A^TA``.
"""

4.4 μs

5.78018

5.78018

5.78018

xxxxxxxxxx
 
opnorm(M,2), opnorm(M),sqrt(maximum(eigvals(M'*M)))

1.1 s

$| | A | |_{\infty} = {max}_{i = 1}^{n} \sum_{j = 1}^{n} | a_{i j} | = | | A^{T} | |_{1}$ maximum of row sums of absolute values of entries

xxxxxxxxxx
 
md"""
- ``||A||_\infty= \max_{i=1}^n \sum_{j=1}^n |a_{ij}| = ||A^T||_1`` maximum of row sums of absolute values of entries
"""

3.4 μs

8.0

8.0

xxxxxxxxxx
 
opnorm(M,Inf),opnorm(M',1)

112 ms

Condition number and error propagation

Solve $A x = b$ , where $b$ is inexact
Let $Δ b$ be the error in $b$ and $Δ x$ be the resulting error in $x$ such that

$A (x + Δ x) = b + Δ b .$

Since $A x = b$ , we get $A Δ x = Δ b$
Therefore $Δ x = A^{- 1} Δ b$

$| | A | | \cdot | | x | | \geq | | b | |$

$| | Δ x | | \leq | | A^{- 1} | | \cdot | | Δ b | |$

$\Rightarrow \frac{| | Δ x | |}{| | x | |} \leq κ (A) \frac{| | Δ b | |}{| | b | |}$

where $κ (A) = | | A | | \cdot | | A^{- 1} | |$ is the condition number of $A$ .

This means that the relative error in the solution is proportional to the relative error of the right hand side. The proportionality factor $κ (A)$ is usually larger (and in most relevant cases significantly larger) than one. Just remark that this estimates does not assume inexact arithmetics.

xxxxxxxxxx
 
md"""
## Condition number and error propagation
​
- Solve ``Ax=b``, where ``b`` is inexact
- Let ``\Delta b`` be the error in ``b`` and ``\Delta x`` be    the resulting error in ``x`` such that
   
     ``A(x+\Delta x)=b+\Delta b.`` 
 
- Since ``Ax=b``, we get ``A\Delta x=\Delta b``
- Therefore   ``\Delta x=A^{-1} \Delta b``
       
   ``||A||\cdot ||x||\geq||b||`` 
​
   ``||\Delta x||\leq ||A^{-1}||\cdot ||\Delta b||``
​
   ``\Rightarrow
   \frac{||\Delta x||}{||x||}
   \leq \kappa(A) \frac{||\Delta b||}{||b||}``
 
 where ``\kappa(A)= ||A||\cdot ||A^{-1}||`` is the *condition number* of ``A``.
​
This means that the relative error in the solution is proportional to the relative error of the right hand side. The proportionality factor ``\kappa(A)`` is usually larger (and in most relevant cases significantly larger) than one. Just remark that
this estimates does not assume inexact arithmetics.
"""

1.5 ms

Let us have an example. We use rational arithmetics in order to perform exact calulations.

xxxxxxxxxx
 
md"""
Let us have an example. We use rational arithmetics in order to perform
exact calulations.
"""

5.3 μs

T_test

Float64

xxxxxxxxxx
 
T_test=Float64

1.0 μs

1.0e6

xxxxxxxxxx
 
a=T_test(1_000_000)

5.5 μs

pert_b

1.0e-9

xxxxxxxxxx
 
pert_b=T_test(1//1_000_000_000)

16.9 ms

2×2 Array{Float64,2}:
 1.0    -1.0
 1.0e6   1.0e6

xxxxxxxxxx
 
A=[ 1 -1;
    a  a]

37.8 ms

1.0e6

xxxxxxxxxx
 
κ=opnorm(A)*opnorm(inv(A))

174 ms

2×2 Array{Float64,2}:
  0.5  5.0e-7
 -0.5  5.0e-7

xxxxxxxxxx
 
inv(A)

9.6 μs

Assume a solution vector:

xxxxxxxxxx
 
md"""
Assume a solution vector:
"""

2.9 μs

Int641

xxxxxxxxxx
 
x=[1,1]

2.4 μs

Create corresponding right hand side:

xxxxxxxxxx
 
md"""
Create corresponding right hand side:
"""

2.7 μs

Float641

0.0

2.0e6

xxxxxxxxxx
 
b=A*x

80.4 ms

Define a perturbation of the right hand side:

xxxxxxxxxx
 
md"""
Define a perturbation of the right hand side:
"""

2.5 μs

Δb

Float641

1.0e-9

1.0e-9

xxxxxxxxxx
 
Δb=[pert_b, pert_b]

3.8 μs

Calculate the error with respect to the unperturbed solution:

xxxxxxxxxx
 
md"""
Calculate the error with respect to the unperturbed solution:
"""

2.4 μs

Δx

Float641

5.0e-10

-5.0e-10

xxxxxxxxxx
 
Δx=inv(A)*(b+Δb)-x

116 ms

Relative error of right hand side:

xxxxxxxxxx
 
md"""
Relative error of right hand side:
"""

2.3 μs

δb

7.071067811865475e-16

xxxxxxxxxx
 
δb=norm(Δb)/norm(b)

7.2 μs

Relative error of solution:

xxxxxxxxxx
 
md"""
Relative error of solution:
"""

2.2 μs

δx

5.000000413703827e-10

xxxxxxxxxx
 
δx=norm(Δx)/norm(x)

69.1 ms

Comparison with condition number based estimate:

xxxxxxxxxx
 
md"""
Comparison with condition number based estimate:
"""

2.3 μs

1.0e6

7.07107e5

xxxxxxxxxx
 
κ,δx/δb

3.2 μs

Complexity: "big O notation"

Let $f, g : V \to R^{+}$ be some functions, where $V = N$ or $V = R$ .

Write

$f (x) = O (g (x)) (x \to \infty)$

if there exists a constant $C > 0$ and $x_{0} \in V$ such that $\forall x > x_{0}, | f (x) | \leq C | g (x) |$

Often, one skips the part " $(x \to \infty)$ "

Examples:

Addition of two vectors: $O (N)$
Matrix-vector multiplication (for matrix where all entries are assumed to be nonzero): $O (N^{2})$

xxxxxxxxxx
 
md"""
## Complexity: "big O notation"
​
  Let ``f,g: \mathbb V \to \mathbb R^+`` be some functions, where ``\mathbb V=\mathbb N`` or ``\mathbb V=\mathbb R``.
​
  Write 
​
```math
  f(x)=O(g(x)) \quad (x\to\infty) 
```
if there exists a constant ``C>0`` and ``x_0\in \mathbb V`` such that
     `` \forall x>x_0, \quad |f(x)|\leq C|g(x)|`` 
​
 Often, one skips the part "``(x \to \infty)``"
 
Examples:
  - Addition of two vectors: ``O(N)``
  - Matrix-vector multiplication (for matrix where all entries are assumed to be nonzero): ``O(N^2)``
​
"""

9.7 μs

A Direct method: Cramer's rule

Solve $A x = b$ by Cramer's rule:

$x_{i} = | \begin{matrix} a_{11} & a_{12} & \dots & a_{1 i - 1} & b_{1} & a_{1 i + 1} & \dots & a_{1 N} \\ a_{21} & \dots & b_{2} & \dots & a_{2 N} \\ ⋮ & ⋮ & ⋮ \\ a_{N 1} & \dots & b_{N} & \dots & a_{N N} \end{matrix} | / | A | (i = 1 \dots N)$

This takes $O (N!)$ operations...

xxxxxxxxxx
 
md"""
## A Direct method: Cramer's rule
​
Solve $Ax=b$ by  Cramer's rule:
​
```math
x_i=\left|
    \begin{matrix}
      a_{11}&a_{12}&\ldots&a_{1i-1}&b_1&a_{1i+1}&\ldots&a_{1N}\\
      a_{21}&      &\ldots&        &b_2&        &\ldots&a_{2N}\\
      \vdots&      &      &        &\vdots&     &      &\vdots\\
      a_{N1}&      &\ldots&        &b_N&        &\ldots&a_{NN}
    \end{matrix}\right|
    / |A| \quad (i=1\dots N)
```
​
This takes     $O(N!)$ operations...
​
"""

5.2 μs

LU decomposition

Gaussian elimination

So let us have a look at Gaussian elimination to solve $A x = b$ . The elementary matrix manipulation step in Gaussian elimination ist the multiplication of row k by $- a_{j k} / a_{k k}$ and its addition to row j such that element $_{j k}$ in the resulting matrix becomes zero. If this is done at once for all $j > k$ , we can express this operation as the left multiplication of $A$ by a lower triangular Gauss transformation matrix $M (A, k)$ .

xxxxxxxxxx
 
md"""
## LU decomposition
​
#### Gaussian elimination 
So let us have a look at Gaussian elimination to solve ``Ax=b``.
The elementary matrix manipulation step in Gaussian elimination ist the multiplication of row k by ``-a_{jk}/a_{kk}`` and its addition to row j such that element
``_{jk}`` in the resulting matrix becomes zero. If this is done at once for all ``j>k``, we can express this operation as the left multiplication of ``A`` by a lower triangular Gauss transformation matrix ``M(A,k)``.  
"""

4.1 μs

xxxxxxxxxx
 
function gausstransform(A,k)
    n=size(A,1)
    M=Matrix(one(eltype(A))I,n,n)
    for j=k+1:n
        M[j,k]=-A[j,k]/A[k,k]
    end
    M
end;

27.6 μs

Define the number type for the following examples:

xxxxxxxxxx
 
md"""
Define the number type for the following examples:
"""

2.1 μs

T_lu

Rational

xxxxxxxxxx
 
T_lu=Rational

2.8 μs

Define a test matrix:

xxxxxxxxxx
 
md"""
Define a test matrix:
"""

2.2 μs

4×4 Array{Rational,2}:
 2//1  1//1  3//1  4//1
 5//1  6//1  7//1  8//1
 7//1  6//1  8//1  5//1
 3//1  4//1  2//1  2//1

xxxxxxxxxx
 
A1=T_lu[2 1 3 4;
        5 6 7 8;
        7 6 8 5;
        3 4 2 2;]   

24.8 ms

This is the Gauss transform for first column:

xxxxxxxxxx
 
md"""
This is the Gauss transform for  first column:
"""

2.3 μs

4×4 Array{Rational{Int64},2}:
  1//1  0//1  0//1  0//1
 -5//2  1//1  0//1  0//1
 -7//2  0//1  1//1  0//1
 -3//2  0//1  0//1  1//1

xxxxxxxxxx
 
gausstransform(A1,1)

71.1 ms

Applying it then sets all elements in the first column to zero besides of the main diagonal element:

xxxxxxxxxx
 
md"""
Applying it then sets all elements in the first column to zero besides of the main diagonal element:
"""

2.2 μs

4×4 Array{Rational,2}:
 2//1  1//1   3//1   4//1
 0//1  7//2  -1//2  -2//1
 0//1  5//2  -5//2  -9//1
 0//1  5//2  -5//2  -4//1

xxxxxxxxxx
 
U1=gausstransform(A1,1)*A1

328 ms

We can repeat this with the second column:

xxxxxxxxxx
 
md"""
We can repeat this with the second column:
"""

2.3 μs

4×4 Array{Rational,2}:
 2//1  1//1    3//1    4//1
 0//1  7//2   -1//2   -2//1
 0//1  0//1  -15//7  -53//7
 0//1  0//1  -15//7  -18//7

xxxxxxxxxx
 
U2=gausstransform(U1,2)*U1

46.7 μs

And the third column:

xxxxxxxxxx
 
md"""
And the third column:
"""

2.2 μs

4×4 Array{Rational,2}:
 2//1  1//1    3//1    4//1
 0//1  7//2   -1//2   -2//1
 0//1  0//1  -15//7  -53//7
 0//1  0//1    0//1    5//1

xxxxxxxxxx
 
U3=gausstransform(U2,3)*U2

34.7 μs

And here, we arrived at a triangular matrix. In the standard Gaussian elimination we would have manipulated the right hand side accordingly.

From here on we would start the backsubstitution which in fact is the solution of a triangular system of equations.

However, let's have a look at what we have done here: we arrived at

$\begin{aligned} U_{1} & = M (A, 1) A \\ U_{2} & = M (U_{1}, 2) U_{1} = M (U_{1}, 2) M (A, 1) A \\ U_{3} & = M (U_{2}, 3) U_{2} = M (U_{2}, 3) M (U_{1}, 2) M (A, 1) A \end{aligned}$

Thus, $A = L U$ with $L = M (A, 1)^{- 1} M (U_{1}, 2)^{- 1} M (U_{2}, 3)^{- 1}$ and $U = U_{3}$ $L$ is a lower triangular matrix and $U$ is a upper triangular matrix.

xxxxxxxxxx
 
md"""
And here, we arrived at a triangular matrix. In the standard Gaussian elimination we would have manipulated the right hand side accordingly. 
​
From here on we would start the backsubstitution which in fact is the solution of a triangular system of equations.
​
However, let's have a look at what we have done here: we arrived at
```math
\begin{align}
U_1&=M(A,1)A\\
U_2&=M(U_1,2)U_1=M(U_1,2)M(A,1)A\\
U_3&=M(U_2,3)U_2=M(U_2,3)M(U_1,2)M(A,1)A
\end{align}
```
Thus, ``A=LU`` with ``L=M(A,1)^{-1}M(U_1,2)^{-1}M(U_2,3)^{-1}`` and ``U=U_3``
``L`` is a lower triangular matrix and ``U`` is a upper triangular matrix.
​
"""

3.9 μs

A first LU decomposition

We can put this together into a function:

xxxxxxxxxx
 
md"""
#### A first LU decomposition
​
We can put this together into a function:
"""

3.0 μs

xxxxxxxxxx
 
function my_first_lu_decomposition(A)
    n=size(A,1)
    L=Matrix(one(eltype(A))I,n,n) # L=I
    U=A
    for k=1:n-1
       M=gausstransform(U,k)
       L=L*inv(M)
       U=M*U
    end
    L,U
end;

31.6 μs

4×4 Array{Rational{Int64},2}:
 1//1  0//1  0//1  0//1
 5//2  1//1  0//1  0//1
 7//2  5//7  1//1  0//1
 3//2  5//7  1//1  1//1

4×4 Array{Rational,2}:
 2//1  1//1    3//1    4//1
 0//1  7//2   -1//2   -2//1
 0//1  0//1  -15//7  -53//7
 0//1  0//1    0//1    5//1

xxxxxxxxxx
 
Lx,Ux=my_first_lu_decomposition(A1)

918 ms

Check for correctness:

xxxxxxxxxx
 
md"""
Check for correctness:
"""

2.4 μs

4×4 Array{Rational{Int64},2}:
 0//1  0//1  0//1  0//1
 0//1  0//1  0//1  0//1
 0//1  0//1  0//1  0//1
 0//1  0//1  0//1  0//1

xxxxxxxxxx
 
Lx*Ux-A1

103 ms

So now we can write $A = L U$ . Solving $L U x = b$ then amounts to solve two triangular systems:

$\begin{aligned} L y & = b \\ U x & = y \end{aligned}$

xxxxxxxxxx
 
md"""
So now we can write  ``A=LU``. 
Solving ``LUx=b`` then amounts to solve two triangular systems:
```math
\begin{align}
   Ly&=b\\
   Ux&=y
\end{align}
```
"""

2.6 μs

my_first_lu_solve (generic function with 1 method)

xxxxxxxxxx
 
function my_first_lu_solve(L,U,b)
   y=inv(L)*b
   x=inv(U)*y
end

17.9 μs

Int641

xxxxxxxxxx
 
b1=[1,2,3,4]

2.3 μs

Rational{Int64}1

182//75

-7//75

-154//75

3//5

xxxxxxxxxx
 
x1=my_first_lu_solve(Lx,Ux,b1)

204 ms

Check...

xxxxxxxxxx
 
md"""
Check...
"""

2.2 μs

Rational{Int64}1

0//1

0//1

0//1

0//1

xxxxxxxxxx
 
A1*x1-b1

147 ms

... in order to be didactical, in this example, we made a very inefficient implementation by creating matrices in each step. We even cheated by using inv in order to solve a triangular system.

xxxxxxxxxx
 
md"""
... in order to be didactical, in this example, we made a very inefficient implementation by creating matrices in each step. We even cheated by using `inv` in order to solve a triangular system.
"""

2.9 μs

A reasonable implementation

Doolittles method (Adapted from wikipedia: LU_decomposition)
This allows to perfrom LU decomposition in-place.

xxxxxxxxxx
 
md"""
#### A reasonable implementation
​
- Doolittles method  (Adapted from [__wikipedia: LU_decomposition__](https://en.wikipedia.org/wiki/LU_decomposition#MATLAB_code_examples))
- This allows to perfrom LU decomposition in-place.
"""

6.8 μs

better_lu_decomposition! (generic function with 1 method)

xxxxxxxxxx
 
function better_lu_decomposition!(LU)
    n = size(LU,1)
    # decomposition of matrix, Doolittle’s Method
    for i = 1:n
        for j = 1:(i - 1)
            alpha = LU[i,j];
            for k = 1:(j - 1)
                alpha = alpha - LU[i,k]*LU[k,j];
            end
            LU[i,j] = alpha/LU[j,j];
        end
        for j = i:n
            alpha = LU[i,j];
            for k = 1:(i - 1)
                alpha = alpha - LU[i,k]*LU[k,j];
            end
            LU[i,j] = alpha;
        end
    end
​
end

42.7 μs

better_lu_solve (generic function with 1 method)

xxxxxxxxxx
 
function better_lu_solve(LU,b)
    n = length(b);
    x = zeros(eltype(LU),n);
    y = zeros(eltype(LU),n);
    # LU= L+U-I
    # find solution of Ly = b
    for i = 1:n
        alpha = zero(eltype(LU));
        for k = 1:i
            alpha = alpha + LU[i,k]*y[k];
        end
        y[i] = b[i] - alpha;
    end
    # find solution of Ux = y
    for i = n:-1:1
        alpha = zero(eltype(LU));
        for k = (i + 1):n
            alpha = alpha + LU[i,k]*x[k];
        end
        x[i] = (y[i] - alpha)/LU[i, i];
    end    
    x
end

40.0 μs

We can then implement a method for linear system solution:

xxxxxxxxxx
 
md"""
We can then implement a method for linear system solution:
"""

2.3 μs

better_solve (generic function with 1 method)

xxxxxxxxxx
 
function better_solve(A,b)
    LU=copy(A)
    better_lu_decomposition!(LU)
    better_lu_solve(LU,b)
end

14.8 μs

Rational1

182//75

-7//75

-154//75

3//5

xxxxxxxxxx
 
x2=better_solve(A1,b1)

88.9 ms

Rational{Int64}1

0//1

0//1

0//1

0//1

xxxxxxxxxx
 
A1*x2-b1

80.6 ms

Pivoting

So far, we ignored the possibility that a diagonal element becomes zero during the LU factorization procedure.

Pivoting tries to remedy the problem that during the algorithm, diagonal elements can become zero. Before undertaking the next Gauss transformation step, we can exchange rows such that we always dividy by the largest of the remaining diagonal elements. This would then in fact result in a decompositon

$P A = L U$

where $P$ is a permutation matrix which can be stored in an integer vector. This approach is called "partial pivoting". Full pivoting in addition would perform column permutations. This would result in another permutation matrix $Q$ and the decomposition

$P A Q = L U$

Almost all practically used LU decomposition implementations use partial pivoting.

xxxxxxxxxx
 
md"""
### Pivoting
So far, we ignored the possibility that a diagonal element becomes zero during the LU factorization procedure.
​
Pivoting tries to remedy the problem that  during the algorithm, diagonal elements can become zero. Before undertaking the next Gauss transformation step, we can exchange rows such that we always dividy by the largest of the remaining diagonal elements.
This would then in fact result in a decompositon
```math
PA=LU
```
where ``P`` is a permutation matrix which can be stored in an integer vector.  This approach is called "partial pivoting". Full pivoting in addition would perform column
permutations. This would result in another permutation matrix ``Q`` and the decomposition
```math
PAQ=LU
```
Almost all practically used LU decomposition implementations use partial pivoting.
"""

4.5 μs

LU Factorization from Julia library

Julia implements a pivoting LU factorization

xxxxxxxxxx
 
md"""
### LU Factorization from Julia library
​
Julia implements a pivoting LU factorization
"""

4.7 μs

lu1

LU{Rational,Array{Rational,2}}
L factor:
4×4 Array{Rational,2}:
 1//1   0//1    0//1  0//1
 5//7   1//1    0//1  0//1
 3//7   5//6    1//1  0//1
 2//7  -5//12  -1//2  1//1
U factor:
4×4 Array{Rational,2}:
 7//1   6//1   8//1    5//1
 0//1  12//7   9//7   31//7
 0//1   0//1  -5//2  -23//6
 0//1   0//1   0//1    5//2

xxxxxxxxxx
 
lu1=lu(A1)

76.3 ms

Like in matlab, the backslash opertor "solves", in this case it solves the LU factorization:

xxxxxxxxxx
 
md"""
Like in matlab, the backslash opertor "solves", in this case it solves the LU factorization:
"""

2.4 μs

Rational{Int64}1

182//75

-7//75

-154//75

3//5

xxxxxxxxxx
 
lu1\b1

114 ms

Of course we can apply \ directly to a matrix. However, behind this always LU decomposition and LU solve are called:

xxxxxxxxxx
 
md"""
Of course we can apply `\` directly to a matrix. However, behind this always LU decomposition and LU solve are called:
"""

2.5 μs

Rational{Int64}1

182//75

-7//75

-154//75

3//5

xxxxxxxxxx
 
x3=A1\b1

792 ms

Rational{Int64}1

0//1

0//1

0//1

0//1

xxxxxxxxxx
 
A1*x3-b1

13.9 μs

LU vs. inv

xxxxxxxxxx
 
md"""
### LU vs. inv
"""

2.3 μs

In principle we could work with the inverse matrix as well:

xxxxxxxxxx
 
md"""
In principle we could work with the inverse matrix as well:
"""

2.3 μs

A1inv

4×4 Array{Rational{Int64},2}:
  68//75  -2//3   2//25   49//75
 -43//75   1//3  -2//25    1//75
 -46//75   1//3   6//25  -53//75
   2//5    0//1  -1//5     1//5

xxxxxxxxxx
 
A1inv=inv(A1)

25.2 μs

Rational{Int64}1

182//75

-7//75

-154//75

3//5

xxxxxxxxxx
 
A1inv*b1

2.6 ms

However, inversion is more complex than the LU factorization.

xxxxxxxxxx
 
md"""
However, inversion is more complex than the LU factorization.
"""

2.5 μs

Some performance tests.

We generate matrices

xxxxxxxxxx
 
md"""
### Some performance tests. 
​
We generate matrices 
"""

3.3 μs

rand_Ab (generic function with 1 method)

xxxxxxxxxx
 
rand_Ab(n)=(100.0I(n)+rand(n,n),rand(n))

18.2 μs

xxxxxxxxxx
 
function perftest_lu(n)
    A,b=rand_Ab(n)
    @elapsed A\b
end;

62.9 μs

perftest_inv (generic function with 1 method)

xxxxxxxxxx
 
function perftest_inv(n)
    A,b=rand_Ab(n)
    @elapsed inv(A)*b
end

55.6 μs

perftest_better (generic function with 1 method)

xxxxxxxxxx
 
function perftest_better(n)
    A,b=rand_Ab(n)
    @elapsed better_solve(A,b)
end

23.2 μs

test_and_plot (generic function with 1 method)

xxxxxxxxxx
 
function test_and_plot(pmax)
    N= 2 .^collect(5:pmax)
    t_inv=perftest_inv.(N)
    t_lu=perftest_lu.(N)
​
    clf()
    loglog(N,t_inv,"-o",label="inv(A)*b")
    loglog(N,t_lu,"-o",label="A\\b")
    loglog(N,1.0e-9*N.^2.75,"k--",label="O(\$N^{2.75}\$)")
    if pmax<12
      t_b=perftest_better.(N)
       loglog(N,t_b,"-o",label="\"better\"")
       loglog(N,1.0e-9*N.^3,"k-",label="O(\$N^{3}\$)")
    end
    xlabel("Number of unknowns N")
    ylabel("Execution time t/s")
    title("Experimental complexity of dense linear system solution")
    grid()
    legend(loc="upper left")
    gcf()
end

44.7 μs

xxxxxxxxxx
 
test_and_plot(11)

7.5 s

The overall complexity in this experiment is around $O (N^{2.75})$ which is in the region of some theoretical estimates.
A good implementation is hard to get right, straightforward code performs worse than the system implementation
Using inversion instead of \ is significantly slower (log scale in the plot!)
For standard floating point types, Julia uses highly optimized versions of LINPACK and BLAS
- Same for python/numpy and many other coding environments

8.3 μs