Notebook 19 – Math 2121, Fall 2020

Below are functions and examples covering the most important algorithms from the whole course.

8.1 μs

xxxxxxxxxx
 
using LinearAlgebra

486 μs

DEFAULT_TOLERANCE

1.0e-9

xxxxxxxxxx
 
DEFAULT_TOLERANCE = 1e-9

40.0 ns

rowop_swap (generic function with 1 method)

xxxxxxxxxx
 
begin
    # the functions here are all helper methods used internally 
    # for the implementation of RREF, pivot_positions, etc.
    
    # it's not recommended to use these functions directly
​
    function find_leading(A, row=1, tol=DEFAULT_TOLERANCE)
        (m, n) = size(A)
        for col=1:n
            if abs(A[row, col]) > tol
                return col
            end
        end
        return 0
    end
    
    function find_nonzeros_in_column(A, row_start, col, tol=DEFAULT_TOLERANCE)
        (m, n) = size(A)
        return [i for i=row_start:m if abs(A[i, col]) > tol]
    end
    
    function rowop_replace(A, source_row, target_row, scalar_factor)
        @assert source_row != target_row
        for j=1:size(A)[2]
            A[target_row,j] += A[source_row,j] * scalar_factor
        end
        return A
    end
    
    function rowop_scale(A, target_row, scalar_factor)
        @assert scalar_factor != 0
        for j=1:size(A)[2]
            A[target_row,j] *= scalar_factor
        end
        A
    end
    
    function rowop_swap(A, s, t)
        for j=1:size(A)[2]
            A[t,j], A[s,j] = A[s,j], A[t,j]
        end
        A
    end
end

111 μs

RREF (generic function with 2 methods)

xxxxxxxxxx
 
function RREF(A, tol=DEFAULT_TOLERANCE)
    # returns reduced echelon form of input matrix A
    
    A = float(copy(A))
    (m, n) = size(A)
​
    row = 1
    for col=1:n
        nonzeros = find_nonzeros_in_column(A, row, col, tol)
        if length(nonzeros) == 0
            continue
        end
​
        i = nonzeros[1]
        for j=nonzeros[2:end]
            if abs(A[i, col]) < abs(A[j, col])
                i = j
            end
        end
​
        if abs(A[i, col] - 1) < tol
            A[i, col] = 1
        else
            A = rowop_scale(A, i, 1 / A[i, col])
        end
​
        for j=nonzeros
            if i != j
                A = rowop_replace(A, i, j, -A[j, col])
            end
        end
​
        if i != row
            A = rowop_swap(A, row, i)
        end
​
        row += 1
    end
​
    for row=m:-1:1
        l = find_leading(A, row, tol)
        if l > 0
            for k=1:row - 1
                if abs(A[k,l]) > tol
                    rowop_replace(A, row, k, -A[k,l])
                end
            end
        end
    end
​
    # round entries to give nicer output; could omit this step
    for row=1:m
        for col=1:n
            if abs(A[row,col] - round(A[row,col])) < tol
                A[row,col] = round(A[row,col])
            end
        end
    end
​
    return A
end

81.6 μs

pivot_positions (generic function with 2 methods)

xxxxxxxxxx
 
function pivot_positions(A, tol=DEFAULT_TOLERANCE)
    # returns list of pivot positions in input matrix A
    rref = RREF(A, tol)
    return [
        (i, find_leading(rref, i)) 
        for i=1:size(A)[1] if find_leading(rref, i) > 0
    ]
end

56.4 μs

solve_linear_system (generic function with 2 methods)

xxxxxxxxxx
 
function solve_linear_system(A, b, tol=DEFAULT_TOLERANCE)
    # returns a solution to Ax = b if one exists, otherwise throws as error
    M = hcat([A b])
    rref = RREF(M)
    pivots = pivot_positions(M, tol)
    n = size(A)[2]
    
    solution_exists = ! ((n + 1) in [j for (i, j) = pivots])
    @assert solution_exists
    
    sol = zeros(n, 1)
    for (i, j) = pivots
        sol[j] = rref[i, end]
    end
    return sol
end

59.1 μs

column_space_basis (generic function with 2 methods)

xxxxxxxxxx
 
function column_space_basis(A, tol=DEFAULT_TOLERANCE)
    # returns a matrix whose columns are a basis for the column space A
    # the columns in the returned matrix are a subset of the columns of A
    pivots = pivot_positions(A, tol)
    ans = zeros(size(A)[1], 0)
    for (i, j) = pivots
        ans = hcat(ans, A[1:end, j])
    end
    return ans
end

42.2 μs

null_space_basis (generic function with 2 methods)

xxxxxxxxxx
 
function null_space_basis(A, tol=DEFAULT_TOLERANCE)
    # returns a matrix whose columns are a basis for the nullspace of A
    
    n = size(A)[2]
    pivots = pivot_positions(A, tol)
    rref = RREF(A, tol)
    free = [j for j=1:n if ! any(j == pivot[2] for pivot=pivots)]
    
    ans = zeros(n, length(free))
    for k = 1:length(free)
        ans[free[k], k] = 1
        for (i, j) = pivots
            ans[j, k] = -rref[i, free[k]]
        end
    end
    return ans
end

72.2 μs

is_diagonalizable (generic function with 1 method)

xxxxxxxxxx
 
function is_diagonalizable(A)
    # returns true if A is diagonalizable
    m, n = size(A)
    return m == n && rank(eigen(A).vectors) == m
end

28.9 μs

gram_schmidt_process (generic function with 1 method)

xxxxxxxxxx
 
function gram_schmidt_process(A; orthonormal=false)
    # returns matrix whose columns are the output of the Gram-Schmidt process
    # applied to the columns of the input matrix `A`. If the parameter `orthonormal` 
    # is set to true, then these columns are normalized to have unit length
​
    m, n = size(A)
    Q = zeros(m, n)
    
    for j=1:n
        a, u = A[1:m, j], A[1:m, j]
        for i = 1:(j - 1)
            v = Q[1:m, i]
            u -= v .* dot(a, v) / dot(v, v)
        end
        Q[1:m, j] = orthonormal ? u / norm(u) : u
    end
​
    return Q
end

90.8 μs

orthogonal_projection (generic function with 1 method)

xxxxxxxxxx
 
function orthogonal_projection(v, A)
    # returns orthogonal projection of v onto Col A
    basis = column_space_basis(A)
    Q = gram_schmidt_process(basis; orthonormal=true)
    return Q * transpose(Q) * v
end

40.7 μs

least_squares_solution (generic function with 2 methods)

xxxxxxxxxx
 
function least_squares_solution(A, b, tol=DEFAULT_TOLERANCE)
    # returns a single least-squares solution to Ax = b
    return solve_linear_system(transpose(A) * A, transpose(A) * b, tol)
end

44.5 μs

pseudoinverse (generic function with 2 methods)

xxxxxxxxxx
 
function pseudoinverse(A, tol=DEFAULT_TOLERANCE)
    U, singular_values, V = svd(A, full=true)
    
    S = transpose(zeros(size(A)))
    for i=1:length(singular_values)
        if abs(singular_values[i]) > tol
            S[i, i] = 1 / singular_values[i]
        end
    end
​
    return V * S * transpose(U)
end         

67.9 μs

random_matrix (generic function with 2 methods)

xxxxxxxxxx
 
function random_matrix(m, n, rank=-1)
    # returns a random m-by-n matrix with the given `rank`. 
    # if the input `rank` is not set, then returned matrix will
    # have full rank which is min(m, n)
    
    # the function is mostly for examples
    
    # there is a tiny probability that the returned 
    # matrix will not have the expected rank
    
    rank = rank == -1 ? minimum([m n]) : rank
    @assert 0 <= rank <= minimum([m n])
    ans = zeros(m, n)
    for i=1:rank
        ans += rand(-1:0.0001:1, m, 1) * rand(-1:0.0001:1, 1, n)
    end
    return 10 * ans
end

47.2 μs

Variable names in Pluto

Pluto notebooks do not allow you to redefine variables names, so a variable can only be assigned once. You won't have this problem if you just import this notebook into Julia without using Pluto, which you can do by opening Julia in your terminal from the directory that contains this file and then running

include("19_Math2121_Fall2020.jl")

To make it easier to extend the examples in this notebook, the variable names used below all begin with underscores, like _A, _b, etc., so as not to crowd the usual namespace of variables. There is no need to continue this convention if you add to this notebook.

9.8 μs

Numerical tolerance

Warning: it is a tricky problem to determine when floating point operations should give the value $0.0$ . Many of the numerical algorithms above depend on a tolerance parameter tol that we use to determine when a floating point x number is zero, by checking if abs(x) < tol. Whether or not you get sensible output from these functions may depend on the value of tol. You can set this manually by passing in a value for tol. For example:

RREF(A) # uses default tolerance

RREF(A, 1e-10) # uses tol = 1e-10

That said, most of the time we don't have to worry about this parameter.

By default tol is set to 1.0e-9.

23.7 μs

Matrix operations

4.0 μs

4×4 Array{Float64,2}:
  0.377141  -6.76723  5.62943    7.26387
  2.51621    1.98108  0.956934  -0.534678
  2.2948    -3.02732  4.62686    4.53796
 -0.785845  -4.35937  2.60611    4.05582

xxxxxxxxxx
 
# creates a 4-by-4 matrix of rank 2
_A = random_matrix(4, 4, 2)

68.1 μs

4×4 Transpose{Float64,Array{Float64,2}}:
  0.377141   2.51621    2.2948   -0.785845
 -6.76723    1.98108   -3.02732  -4.35937
  5.62943    0.956934   4.62686   2.60611
  7.26387   -0.534678   4.53796   4.05582

xxxxxxxxxx
 
transpose(_A)

46.0 ns

xxxxxxxxxx
 
size(_A)

207 ns

xxxxxxxxxx
 
size(_A)[1]

290 ns

xxxxxxxxxx
 
size(_A)[2]

204 ns

4×1 Array{Int64,2}:
  0
  2
  2
 -2

xxxxxxxxxx
 
_b = rand(-2:2, 4, 1)

3.6 μs

4×1 Array{Int64,2}:
 -2
  2
  1
  2

xxxxxxxxxx
 
_c = rand(-2:2, 4, 1)

3.3 μs

4×1 Array{Float64,2}:
 -16.803342999999998
   6.945391
  -5.8768449999999985
 -11.618156999999998

xxxxxxxxxx
 
# matrix multiplication
_A * _b

19.3 μs

4×1 Array{Int64,2}:
 -2
  4
  3
  0

xxxxxxxxxx
 
# vector addition
_b + _c

4.4 μs

4×1 Array{Int64,2}:
  0
  4
  2
 -4

xxxxxxxxxx
 
# multiply corresponding entries of two matrices/vectors of same size together
# note this is NOT matrix multiplication
_b .* _c

101 ms

4×4 Diagonal{Int64,Array{Int64,1}}:
 0  ⋅  ⋅   ⋅
 ⋅  2  ⋅   ⋅
 ⋅  ⋅  2   ⋅
 ⋅  ⋅  ⋅  -2

xxxxxxxxxx
 
# create diagonal matrix from vector
Diagonal(_b[:])

5.0 ms

4×4 Array{Float64,2}:
  0.377141  -6.76723  5.62943    7.26387
  2.51621    1.98108  0.956934  -0.534678
  2.2948    -3.02732  4.62686    4.53796
 -0.785845  -4.35937  2.60611    4.05582

xxxxxxxxxx
 
_A

49.0 ns

2×2 Array{Float64,2}:
 2.51621   1.98108
 2.2948   -3.02732

xxxxxxxxxx
 
# submatrices 
_A[2:3, 1:2] # 2-by-2 submatrix

1.6 μs

Float641

-6.76723

1.98108

-3.02732

-4.35937

xxxxxxxxxx
 
_A[:, 2] # column 2

4.2 μs

Float641

2.51621

1.98108

0.956934

-0.534678

xxxxxxxxxx
 
_A[2, :] # row 2

3.5 μs

-3.0273244999999993

xxxxxxxxxx
 
_A[3, 2] # entry in row 3, column 2

235 ns

_square_matrix

3×3 Array{Int64,2}:
 1  2  3
 3  4  2
 0  1  0

xxxxxxxxxx
 
_square_matrix = [1 2 3; 3 4 2; 0 1 0]

35.2 ms

7.0

xxxxxxxxxx
 
det(_square_matrix)

19.0 μs

xxxxxxxxxx
 
tr(_square_matrix)

233 ns

3×6 Array{Int64,2}:
 1  2  3  1  0  0
 3  4  2  0  1  0
 0  1  0  0  0  1

xxxxxxxxxx
 
hcat(_square_matrix, I)

29.7 μs

3×6 Array{Float64,2}:
 1.0  0.0  0.0  -0.285714   0.428571  -1.14286
 0.0  1.0  0.0   0.0        0.0        1.0
 0.0  0.0  1.0   0.428571  -0.142857  -0.285714

xxxxxxxxxx
 
RREF(hcat(_square_matrix, I))

129 ms

3×3 Array{Float64,2}:
 -0.285714   0.428571  -1.14286
 -0.0        0.0        1.0
  0.428571  -0.142857  -0.285714

xxxxxxxxxx
 
_square_matrix^-1

49.5 μs

3×3 Array{Int64,2}:
  7  13   7
 15  24  17
  3   4   2

xxxxxxxxxx
 
_square_matrix^2

6.4 μs

Row reduction

5.5 μs

4×4 Array{Float64,2}:
  0.377141  -6.76723  5.62943    7.26387
  2.51621    1.98108  0.956934  -0.534678
  2.2948    -3.02732  4.62686    4.53796
 -0.785845  -4.35937  2.60611    4.05582

xxxxxxxxxx
 
_A

40.0 ns

4×4 Array{Float64,2}:
  1.0  0.0   0.991744   0.606026
 -0.0  1.0  -0.776596  -1.03962
  0.0  0.0   0.0        0.0
  0.0  0.0   0.0       -0.0

xxxxxxxxxx
 
RREF(_A)

92.7 ms

Tuple{Int64,Int64}11

xxxxxxxxxx
 
pivot_positions(_A)

37.5 ms

4×2 Array{Float64,2}:
  0.377141  -6.76723
  2.51621    1.98108
  2.2948    -3.02732
 -0.785845  -4.35937

xxxxxxxxxx
 
column_space_basis(_A)

62.4 ms

4×2 Array{Float64,2}:
 -0.991744  -0.606026
  0.776596   1.03962
  1.0        0.0
  0.0        1.0

xxxxxxxxxx
 
null_space_basis(_A)

69.3 ms

4×2 Array{Float64,2}:
 0.0           0.0
 0.0           1.11022e-16
 8.88178e-16   0.0
 4.44089e-16  -8.88178e-16

xxxxxxxxxx
 
_A * null_space_basis(_A) # should be close to zero

18.5 μs

xxxxxxxxxx
 
rank(_A)

81.0 μs

false

xxxxxxxxxx
 
# to detect if some vectors are linearly independent, assemble the vectors
# as the columns of a matrix A and check whether rank(A) == size(A)[2]
​
rank(_A) == size(_A)[2] # columns of matrix not linearly dependent

24.1 μs

Solving linear systems

3.4 μs

4×4 Array{Float64,2}:
  -3.57174   -7.35301   -3.42201   3.81875
   5.69388    6.33517   -1.96969   1.70079
   7.69715    3.05182  -10.2608   10.2693
 -11.0038   -13.5657     1.98354  -1.37461

xxxxxxxxxx
 
_M = random_matrix(4, 4, 2)

12.9 μs

_test_vec

4×1 Array{Float64,2}:
  -5.958647520884568
   5.114630507983956
   2.4275068057374454
 -10.96086766946625

xxxxxxxxxx
 
_test_vec = _M * rand(size(_M)[2], 1)

8.1 μs

_test_sol

4×1 Array{Float64,2}:
 -0.007336740627503224
  0.8139325807742661
  0.0
  0.0

xxxxxxxxxx
 
# solving _M * v == _test_vec
_test_sol = solve_linear_system(_M, _test_vec)

43.4 ms

4×1 Array{Float64,2}:
  -5.95865
   5.11463
   2.42751
 -10.9609

4×1 Array{Float64,2}:
  -5.95865
   5.11463
   2.42751
 -10.9609

xxxxxxxxxx
 
# these vectors should be the same, meaning our solution is a valid solution 
_M * _test_sol, _test_vec

5.0 μs

_dim_nul_M

xxxxxxxxxx
 
_dim_nul_M = size(null_space_basis(_M))[2]

11.8 μs

_other_sol

4×1 Array{Float64,2}:
 0.04221002984635365
 0.7791793775967624
 0.06048291925849458
 0.03362373291973686

xxxxxxxxxx
 
# to generate arbitrary solutions to Mx = b, just add to any given solution
# linear combinations of the basis element for Nul(M)
_other_sol = _test_sol + null_space_basis(_M) * rand(_dim_nul_M, 1)

18.2 μs

4×1 Array{Float64,2}:
  -5.95865
   5.11463
   2.42751
 -10.9609

4×1 Array{Float64,2}:
  -5.95865
   5.11463
   2.42751
 -10.9609

xxxxxxxxxx
 
# these should still be equal
_M * _other_sol, _test_vec

6.3 μs

_vec_not_in_column_space

Float641

-0.684311

-0.80944

1.0

0.0

x
 
# beware: if b is not in Col A then A \ b gives nonsensical output or an error
​
# the vector below is the first element of a basis for the null space of A^T
# the null space of A^T is the orthogonal complement of the column space of A
# hence this vector is *not* contained in the column space of A
_vec_not_in_column_space = null_space_basis(transpose(_A))[:, 1]

17.8 μs

x
 
# method `solve_linear_system` will raise an error if no solution exists
# solve_linear_system(_A, _vec_not_in_column_space)

39.0 ns

Complex numbers

5.9 μs

5 + 6im

x
 
5 + 6 * im

203 ns

x
 
real(5 + 6 * im)

208 ns

x
 
imag(5 + 6 * im)

300 ns

5 - 6im

x
 
conj(5 + 6 * im)

271 ns

9 + 8im

x
 
(5 + 6 * im) + (4 + 2 * im)

222 ns

8 + 34im

xxxxxxxxxx
 
(5 + 6 * im) * (4 + 2 * im)

195 ns

1.6 + 0.7im

xxxxxxxxxx
 
(5 + 6 * im) / (4 + 2 * im)

224 ns

7.810249675906654

xxxxxxxxxx
 
abs(5 + 6 * im)

899 ns

true

x
 
abs(5 + 6 * im) == sqrt((5 + 6 * im) * conj(5 + 6 * im))

962 ns

2 + 3im

xxxxxxxxxx
 
_z = 2 + 3 * im

198 ns

2×2 Array{Int64,2}:
 2  -3
 3   2

xxxxxxxxxx
 
[real(_z) -imag(_z); imag(_z) real(_z)]

28.0 ms

Eigenvalues and eigenvectors

2.4 μs

4×4 Array{Float64,2}:
  6.0118   -5.45765   -5.14928  -0.975021
  2.02876   0.464118   4.78987  -3.90178
 -4.63644  -6.61134   -3.49888  -9.10147
 10.7833    4.37568   -1.54777  13.5531

xxxxxxxxxx
 
# random 4-by-4 matrix of rank 3
_B = random_matrix(4, 4, 3) # diagonalizable

14.8 μs

4×4 Array{Int64,2}:
 1  2  3  4
 0  1  2  3
 0  0  1  2
 0  0  0  1

xxxxxxxxxx
 
_C = [1 2 3 4; 0 1 2 3; 0 0 1 2; 0 0 0 1] # not diagonalizable

34.4 ms

true

xxxxxxxxxx
 
is_diagonalizable(_B)

142 ms

false

xxxxxxxxxx
 
is_diagonalizable(_C)

3.4 ms

Complex{Float64}1

-7.42946e-15+0.0im

0.254658-6.67987im

0.254658+6.67987im

16.0208+0.0im

xxxxxxxxxx
 
# eigenvalus of a matrix
eigen(_B).values

46.1 μs

Complex{Float64}1

-4.90617e-15+0.0im

0.254658-6.67987im

0.254658+6.67987im

16.0208+0.0im

 
eigen(transpose(_B)).values # same eigenvalues as _B, some may be complex numbers

55.2 μs

4×4 Array{Complex{Float64},2}:
 -0.172547+0.0im  -0.0297456+0.418441im  -0.0297456-0.418441im   0.260016+0.0im
 -0.734333+0.0im   -0.671807-0.0im        -0.671807+0.0im       -0.289454+0.0im
  0.495264+0.0im    0.118236+0.491568im    0.118236-0.491568im  -0.359243+0.0im
  0.430925+0.0im   0.0936165-0.329111im   0.0936165+0.329111im   0.848264+0.0im

xxxxxxxxxx
 
# eigenvectors of a matrix
eigen(_B).vectors

46.1 μs

4×4 Array{Complex{Float64},2}:
 -0.435685+0.0im  -0.160302+0.442395im   -0.160302-0.442395im     0.761407+0.0im
  0.215249+0.0im   0.591921-0.0im         0.591921+0.0im        0.00637559+0.0im
  0.711081+0.0im   0.440853+0.198845im    0.440853-0.198845im    -0.246823+0.0im
  0.508144+0.0im   0.437822-0.0513943im   0.437822+0.0513943im    0.599414+0.0im

xxxxxxxxxx
 
eigen(transpose(_B)).vectors # different eigenvectors from B

59.8 μs

Float641

1.0

1.0

1.0

1.0

xxxxxxxxxx
 
eigen(_C).values

22.4 μs

4×4 Array{Float64,2}:
 1.0  -1.0           1.0          -1.0
 0.0   1.11022e-16  -1.11022e-16   1.11022e-16
 0.0   0.0           1.2326e-32   -1.2326e-32
 0.0   0.0           0.0           1.36846e-48

xxxxxxxxxx
 
eigen(_C).vectors

22.5 μs

4×4 Array{Complex{Float64},2}:
 -0.172547+0.0im  -0.0297456+0.418441im  -0.0297456-0.418441im   0.260016+0.0im
 -0.734333+0.0im   -0.671807-0.0im        -0.671807+0.0im       -0.289454+0.0im
  0.495264+0.0im    0.118236+0.491568im    0.118236-0.491568im  -0.359243+0.0im
  0.430925+0.0im   0.0936165-0.329111im   0.0936165+0.329111im   0.848264+0.0im

4×4 Diagonal{Complex{Float64},Array{Complex{Float64},1}}:
 -7.42946e-15+0.0im           ⋅                   ⋅                  ⋅    
              ⋅      0.254658-6.67987im           ⋅                  ⋅    
              ⋅               ⋅          0.254658+6.67987im          ⋅    
              ⋅               ⋅                   ⋅          16.0208+0.0im

xxxxxxxxxx
 
# if matrix is diagonalizable
_P, _D = eigen(_B).vectors, Diagonal(eigen(_B).values)

4.1 ms

2.830578330796744e-14

xxxxxxxxxx
 
norm(_P * _D * _P^-1 - _B) # if this is very small, then B equals P * D * P^-1

1.1 s

_symmetric

4×4 Array{Float64,2}:
 12.0236   -3.42889    -9.78571    9.80823
 -3.42889   0.928236   -1.82147    0.473901
 -9.78571  -1.82147    -6.99775  -10.6492
  9.80823   0.473901  -10.6492    27.1062

xxxxxxxxxx
 
# if matrix is symmetric, the eigenvectors will be orthonormal already 
_symmetric = _B + transpose(_B)

2.6 μs

_eigenvectors

4×4 Array{Float64,2}:
 0.323252   0.234151   0.793089   0.4601
 0.196796   0.918942  -0.341333  -0.0175587
 0.911403  -0.24589   -0.12072   -0.307098
 0.161629  -0.200641  -0.489826   0.832881

xxxxxxxxxx
 
_eigenvectors = eigen(_symmetric).vectors

85.4 μs

4×4 Array{Float64,2}:
  1.0          -7.95206e-18   1.78588e-16   5.29273e-17
 -7.95206e-18   1.0           9.42052e-17  -2.45753e-17
  1.78588e-16   9.42052e-17   1.0          -2.18036e-16
  5.29273e-17  -2.45753e-17  -2.18036e-16   1.0

xxxxxxxxxx
 
# shows that eigenvectors are already orthonormal
transpose(_eigenvectors) * _eigenvectors

14.6 μs

Inner products and orthogonal projections

2.3 μs

3.4641016151377544

xxxxxxxxxx
 
# length of a vector
norm(_b) # returns the same thing as sqrt(sum(_b .* _b))

691 ns

xxxxxxxxxx
 
# dot product
dot(_b, _c) # returns the same thing as sum(_b .* _c)

453 ns

_basis

5×3 Array{Float64,2}:
 -2.43045    2.55144  -6.17487
 -3.3509     9.80676   4.91765
  0.676911  -5.97453  -7.44102
 -6.03631    5.55339  -7.68659
 -0.555015  -2.87016  -7.55904

xxxxxxxxxx
 
# creates 5-by-3 matrix with rank 3, meaning columns are linearly independent
_basis = random_matrix(5, 3, 3)

13.8 μs

xxxxxxxxxx
 
rank(_basis) # should be 3

24.6 μs

_orthogonal_basis

5×3 Array{Float64,2}:
 -2.43045   -0.80474  -2.92988
 -3.3509     5.17955  -0.557834
  0.676911  -5.03979  -0.0171806
 -6.03631   -2.78207   1.61207
 -0.555015  -3.63658  -1.35563

xxxxxxxxxx
 
_orthogonal_basis = gram_schmidt_process(_basis)

9.9 μs

3×3 Array{Float64,2}:
 54.3388        1.29626e-14  -2.6544e-14
  1.29626e-14  73.8394       -7.95571e-15
 -2.6544e-14   -7.95571e-15  13.3321

xxxxxxxxxx
 
# diagonal entries should be positive and off-diagonal entries close to zero
transpose(_orthogonal_basis) * _orthogonal_basis

6.1 μs

_orthonormal_basis

5×3 Array{Float64,2}:
 -0.32971    -0.0936508  -0.802416
 -0.454575    0.602765   -0.152776
  0.0918283  -0.5865     -0.00470532
 -0.818872   -0.32376     0.441503
 -0.0752921  -0.423203   -0.37127

xxxxxxxxxx
 
# setting orthonormal=true will return matrix with unit length orthogonal columns
_orthonormal_basis = gram_schmidt_process(_basis, orthonormal=true)

4.4 ms

3×3 Array{Float64,2}:
  1.0          -6.35764e-17  3.31234e-16
 -6.35764e-17   1.0          2.36696e-16
  3.31234e-16   2.36696e-16  1.0

xxxxxxxxxx
 
# diagonal entries should be 1.0 and off-diagonal entries close to zero
transpose(_orthonormal_basis) * _orthonormal_basis

7.2 μs

5×1 Array{Int64,2}:
  2
 -2
 -1
 -4
 -1

xxxxxxxxxx
 
_v = rand(-5:5, 5, 1)

3.6 μs

_v_hat

5×1 Array{Float64,2}:
  0.9157012881291711
 -0.6344293450681221
 -0.19998600158597624
 -4.355743010034982
  0.3483621248361716

xxxxxxxxxx
 
_v_hat = orthogonal_projection(_v, _basis)

31.3 μs

5×4 Array{Float64,2}:
  1.0   0.0  0.0   2.82911
  0.0   1.0  0.0   1.27133
 -0.0  -0.0  1.0  -0.736532
  0.0   0.0  0.0  -0.0
  0.0   0.0  0.0   0.0

xxxxxxxxxx
 
# if _v_hat is in the span of _basis, then there should be no pivots in last column
RREF(hcat([_basis _v_hat])) 

8.7 μs

3×1 Array{Float64,2}:
 1.3131522849785959e-15
 3.2528969620086696e-15
 5.634650153790826e-15

xxxxxxxxxx
 
# entries should be close to zero as _v - _v_hat is orthogonal to span of _basis
transpose(_basis) * (_v - _v_hat)

7.2 μs

4×5 Array{Float64,2}:
  1.0  0.0   0.991744   0.606026   0.914466
 -0.0  1.0  -0.776596  -1.03962    0.0978552
  0.0  0.0   0.0        0.0        0.0
  0.0  0.0   0.0       -0.0       -0.0

xxxxxxxxxx
 
RREF(hcat([_A orthogonal_projection(_b, _A)]))

25.7 μs

Least-squares solutions

3.2 μs

_least_squares_solution

4×1 Array{Float64,2}:
 0.9144659461447051
 0.09785518382324851
 0.0
 0.0

xxxxxxxxxx
 
_least_squares_solution = least_squares_solution(_A, _b)

5.2 ms

4×1 Array{Float64,2}:
 0.9144659461447051
 0.09785518382324851
 0.0
 0.0

xxxxxxxxxx
 
# returns same thing as least_squares_solution(_A, _b)
solve_linear_system(transpose(_A) * _A, transpose(_A) * _b)

18.7 μs

4×1 Array{Float64,2}:
 11.1937
  6.62625
  5.95538
 -0.105076

4×1 Array{Float64,2}:
 11.1937
  6.62625
  5.95538
 -0.105076

xxxxxxxxxx
 
transpose(_A) * _A * _least_squares_solution, transpose(_A) * _b

9.6 μs

_dim_nul

xxxxxxxxxx
 
_dim_nul = size(null_space_basis(transpose(_A) * _A))[2]

17.8 μs

_another_least_squares_solution

4×1 Array{Float64,2}:
 -0.31484149960282926
  1.493477642011315
  0.7712839971084404
  0.7662893290500663

xxxxxxxxxx
 
# to form an arbitrary least-squares solution to Ax = b, just add
# to a single least-squares solution any linear combination of a 
# basis for the null space of transpose(A) * A
_another_least_squares_solution = _least_squares_solution + null_space_basis(transpose(_A) * _A) * rand(_dim_nul, 1)

20.7 μs

4×1 Array{Float64,2}:
 -0.317326
  2.49485
  1.80228
 -1.14522

4×1 Array{Float64,2}:
 -0.317326
  2.49485
  1.80228
 -1.14522

xxxxxxxxxx
 
_A * _another_least_squares_solution, _A * _least_squares_solution

5.8 μs

Singular value decompositions

2.5 μs

xxxxxxxxxx
 
_m, _n = 4, 5

64.0 ns

4×5 Array{Float64,2}:
 -3.4233    7.59986    0.0224568  -8.46628    5.98055
 -4.28842  -0.502641   5.19409     1.47641    7.49157
  1.62312  -1.38475   -1.15415     1.33975   -2.83553
 -3.4536    0.860067   3.53104    -0.335714   6.03324

xxxxxxxxxx
 
_X = random_matrix(_m, _n, 2)

13.7 μs

SVD{Float64,Float64,Array{Float64,2}}
U factor:
4×4 Array{Float64,2}:
  0.718735   0.661076   -0.209989   0.0479817
  0.479036  -0.664983   -0.512437  -0.256378
 -0.247517   0.0170652  -0.62418    0.74084
  0.438953  -0.347108    0.551101   0.61897
singular values:
4-element Array{Float64,1}:
 16.002967090628864
 10.105428798819345
  1.1687268383051277e-15
  4.133566008353511e-16
Vt factor:
5×5 Array{Float64,2}:
 -0.401955   0.371292   0.271195  -0.365978   0.702202
  0.179619   0.498363  -0.463561  -0.637207  -0.313766
  0.537101  -0.502761   0.344718  -0.564462   0.145963
 -0.65493   -0.578057  -0.327772  -0.339399  -0.119548
  0.297923  -0.163886  -0.696634   0.161919   0.610628

xxxxxxxxxx
 
_U, _singular_values, _V = svd(_X, full=true)

49.5 ms

4×5 Array{Float64,2}:
 16.003   0.0     0.0          0.0          0.0
  0.0    10.1054  0.0          0.0          0.0
  0.0     0.0     1.16873e-15  0.0          0.0
  0.0     0.0     0.0          4.13357e-16  0.0

xxxxxxxxxx
 
begin
    _Sigma = zeros(size(_X))
    for i=1:length(_singular_values)
        _Sigma[i, i] = _singular_values[i]
    end
    _Sigma
end

1.8 μs

5.393785051703609e-15

xxxxxxxxxx
 
# should be very small value
norm(_U * _Sigma * transpose(_V) - _X)

9.2 μs

_X_plus

5×4 Array{Float64,2}:
 -0.00630251  -0.023852    0.00652033  -0.0171951
  0.0492776   -0.0216802  -0.00490116  -0.00693373
 -0.0181451    0.0386225  -0.00497738   0.0233614
 -0.0581218    0.0309759   0.00458449   0.0118486
  0.0110118    0.0416671  -0.0113908    0.0300385

xxxxxxxxxx
 
_X_plus = pseudoinverse(_X)

28.5 ms

8.40622e-15

5.89168e-17

xxxxxxxxxx
 
# should be very small values 
norm(_X * _X_plus * _X - _X), norm(_X_plus * _X * _X_plus - _X_plus)

8.3 μs