# Example of an ill-conditioned system¶

In :
import LinearAlgebra
import Random
Random.seed!(0)

# Shorthand notation
norm = LinearAlgebra.norm
cond = LinearAlgebra.cond

Out:
cond (generic function with 15 methods)
In :
# # Example of an ill-conditioned system

In :
# Parameters of linear system
A = [1 1; 1 (1-1e-12)]
b = [0; 1e-12]

# Use much more accurate BigFloat format
exact_A = [1 1; 1 (1- BigFloat("1e-12"))]
exact_b = [0, BigFloat("1e-12")]

# Relative error on data
relative_error_b = norm(b - exact_b) / norm(b)
relative_error_A = norm(A - exact_A) / norm(A)
println("Relative error on b: $relative_error_b") println("Relative error on A:$relative_error_A")

Relative error on b: 2.011335237074438503729410008359924320584757100186531629388143594799196078353068e-17
Relative error on A: 1.106086006075478567832268914387623212025069211192183133525679964913483566248121e-17

In :
# Exact and approximate solutions
x_exact = [1; -1]
x_approx = A\b

relative_error_x = norm(x_exact - x_approx)/norm(x_approx)
println("Relative error on x: $relative_error_x")  Relative error on x: 2.2121720121379184e-5  In : # Condition number of A cond_A = cond(A)  Out: 4.000660598318207e12 The relative error on the solution is much larger than the relative error on the data. The ratio between the two is of the same order of magnitude as the condition number of the matrix$A$. # LU decomposition¶ In : function lu(A) n = size(A) L = [i == j ? 1.0 : 0.0 for i in 1:n, j in 1:n] U = copy(A) for i in 1:n-1 for r in i+1:n U[i, i] == 0 && error("Pivotal entry is zero!") ratio = U[r, i] / U[i, i] L[r, i] = ratio U[r, i:end] -= U[i, i:end] * ratio end end return L, U end ns = 2 .^ collect(5:9) times = zeros(length(ns)) for (i, n) in enumerate(ns) A = randn(n, n) result = @timed lu(A) L, U = result.value times[i] = result.time error = norm(L*U - A) println("Error:$(LinearAlgebra.norm(L*U - A)), time: $(result.time)") end  Error: 1.8246083669574607e-12, time: 0.456947637 Error: 9.975383169794571e-13, time: 0.005081298 Error: 3.4400222673534816e-11, time: 0.100179736 Error: 1.2980101641877466e-10, time: 0.737644245 Error: 1.9962768190723483e-9, time: 2.961592128  # LU decomposition with pivoting¶ In : function lu_pivot(A) n = size(A) L, U = zeros(n, 0), copy(A) P = [i == j ? 1.0 : 0.0 for i in 1:n, j in 1:n] function swap_rows!(i, j, matrices...) for M in matrices M_row_i = M[i, :] M[i, :] = M[j, :] M[j, :] = M_row_i end end for i in 1:n-1 # Pivoting index_row_pivot = i - 1 + argmax(abs.(U[i:end, i])) swap_rows!(i, index_row_pivot, U, L, P) # Usual Gaussian transformation c = [zeros(i-1); 1.0; zeros(n-i)] for r in i+1:n ratio = U[r, i] / U[i, i] c[r] = ratio U[r, i:end] -= U[i, i:end] * ratio end L = [L c] end L = [L [zeros(n-1); 1.0]] return P, L, U end ns = 2 .^ collect(5:9) times = zeros(length(ns)) for (i, n) in enumerate(ns) A = randn(n, n) L, U = lu(A) P, Lpivot, Upivot = lu_pivot(A) println("\n-- n =$n --")
println("Condition numbers without pivoting: κ(L) = $(cond(L)), κ(U) =$(cond(U))")
println("Condition numbers with    pivoting: κ(L) = $(cond(Lpivot)), κ(U) =$(cond(Upivot))")
end

-- n = 32 --
Condition numbers without pivoting: κ(L) = 178637.01929136374, κ(U) = 161595.53057380614
Condition numbers with    pivoting: κ(L) = 28.691391551934966, κ(U) = 193.19052697115316

-- n = 64 --
Condition numbers without pivoting: κ(L) = 243901.7856336236, κ(U) = 206804.18071457956
Condition numbers with    pivoting: κ(L) = 83.9559265157693, κ(U) = 402.1768188798624

-- n = 128 --
Condition numbers without pivoting: κ(L) = 5.856999082920966e6, κ(U) = 6.169501793769366e6
Condition numbers with    pivoting: κ(L) = 207.60793737470857, κ(U) = 599.3815911475028

-- n = 256 --
Condition numbers without pivoting: κ(L) = 4.2755305949334687e8, κ(U) = 4.729207635483682e8
Condition numbers with    pivoting: κ(L) = 447.0937209341298, κ(U) = 1009.1733665132496

-- n = 512 --
Condition numbers without pivoting: κ(L) = 1.4221750935427153e8, κ(U) = 1.3737670790079305e8
Condition numbers with    pivoting: κ(L) = 989.5470052618462, κ(U) = 3785.056520009892


# Reducing the bandwidth via permutations¶

In :
import SparseArrays
A = LinearAlgebra.diagm(0 => 10*ones(8), 1 => ones(7), -1 => ones(7), 7 => ones(1), -7 => ones(1))
A = SparseArrays.SparseMatrixCSC(A)
display(A)

# Construct permutation matrix obtained by Cuthill-McKee
σ = [1, 2, 4, 6, 8, 7, 5, 3]
Pσ = SparseArrays.spzeros(8, 8)
for j in 1:8
Pσ[σ[j], j] = 1
end

# Matrix after permutation
new_A = Pσ*A*Pσ'
display(new_A)

8×8 SparseArrays.SparseMatrixCSC{Float64,Int64} with 24 stored entries:
[1, 1]  =  10.0
[2, 1]  =  1.0
[8, 1]  =  1.0
[1, 2]  =  1.0
[2, 2]  =  10.0
[3, 2]  =  1.0
[2, 3]  =  1.0
[3, 3]  =  10.0
[4, 3]  =  1.0
[3, 4]  =  1.0
[4, 4]  =  10.0
[5, 4]  =  1.0
[4, 5]  =  1.0
[5, 5]  =  10.0
[6, 5]  =  1.0
[5, 6]  =  1.0
[6, 6]  =  10.0
[7, 6]  =  1.0
[6, 7]  =  1.0
[7, 7]  =  10.0
[8, 7]  =  1.0
[1, 8]  =  1.0
[7, 8]  =  1.0
[8, 8]  =  10.0
8×8 SparseArrays.SparseMatrixCSC{Float64,Int64} with 24 stored entries:
[1, 1]  =  10.0
[2, 1]  =  1.0
[3, 1]  =  1.0
[1, 2]  =  1.0
[2, 2]  =  10.0
[4, 2]  =  1.0
[1, 3]  =  1.0
[3, 3]  =  10.0
[5, 3]  =  1.0
[2, 4]  =  1.0
[4, 4]  =  10.0
[6, 4]  =  1.0
[3, 5]  =  1.0
[5, 5]  =  10.0
[7, 5]  =  1.0
[4, 6]  =  1.0
[6, 6]  =  10.0
[8, 6]  =  1.0
[5, 7]  =  1.0
[7, 7]  =  10.0
[8, 7]  =  1.0
[6, 8]  =  1.0
[7, 8]  =  1.0
[8, 8]  =  10.0