
import numpy as np 
from numpy.linalg import norm


vector_row = np.array([[1, -5, 3, 2, 4]])
print(vector_row)
new_vector = vector_row.T
print(new_vector)

norm_1 = norm(new_vector, 1)
norm_2 = norm(new_vector, 2)
norm_inf = norm(new_vector, np.inf)

print('L_1 is: %.1f'%norm_1)
print('L_2 is: %.1f'%norm_2)

[[ 1 -5  3  2  4]]
[[ 1]
 [-5]
 [ 3]
 [ 2]
 [ 4]]
L_1 is: 15.0
L_2 is: 7.4


u = np.array([[0, 3, 2]])
v = np.array([[4, 1, 1]])
w = np.array([[0, -2, 0]])

x = 3*u - 2*v + 4*w

print(x)

[[-8 -1  4]]


print(u)
print(v)
print(v.T)

u.dot(v)   #  .... Read the error!

[[0 3 2]]
[[4 1 1]]
[[4]
 [1]
 [1]]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-33f082660c40> in <module>
      3 print(v.T)
      4 
----> 5 u.dot(v)   #  .... Read the error!

ValueError: shapes (1,3) and (1,3) not aligned: 3 (dim 1) != 1 (dim 0)


print(u.shape)
print(v.shape)

(1, 3)
(1, 3)


u = np.array([[0, 3, 2]])
v = np.array([4, 1, 1])    # notice the change here! 

x = u * v 

print(x)

[[0 3 2]]


print(u)
print(v)
print(v.T) # Does nothing ... 

u.dot(v)   # this works here!

[[0 3 2]]
[4 1 1]
[4 1 1]

array([5])


print(u.shape)
print(v.shape)

(1, 3)
(3,)


P = np.array([[1, 7], [2, 3], [5, 0]])
Q = np.array([[2, 6, 3, 1], [1, 2, 3, 4]])
print(P)
print(Q)
print(np.dot(P, Q))

np.dot(Q, P) # can't do that

[[1 7]
 [2 3]
 [5 0]]
[[2 6 3 1]
 [1 2 3 4]]
[[ 9 20 24 29]
 [ 7 18 15 14]
 [10 30 15  5]]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-9-8884e7614908> in <module>
      5 print(np.dot(P, Q))
      6 
----> 7 np.dot(Q, P) # can't do that

<__array_function__ internals> in dot(*args, **kwargs)

ValueError: shapes (2,4) and (3,2) not aligned: 4 (dim 1) != 3 (dim 0)


P.shape, Q.shape

((3, 2), (2, 4))


from numpy.linalg import det

M = np.array(
             [  [0, 2, 1, 3], 
                [3, 2, 8, 1], 
                [1, 0, 0, 3],
                [0, 3, 2, 1]  ]
            )
print('M:\n', M)

print('Determinant: %.1f' % det(M))

I = np.eye(4)  # eye-dentity matrix 

print('I:\n', I)
print('M*I:\n', np.dot(M, I))

M:
 [[0 2 1 3]
 [3 2 8 1]
 [1 0 0 3]
 [0 3 2 1]]
Determinant: -38.0
I:
 [[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]
M*I:
 [[0. 2. 1. 3.]
 [3. 2. 8. 1.]
 [1. 0. 0. 3.]
 [0. 3. 2. 1.]]


from numpy.linalg import inv

print('Inv M:\n', inv(M))

Inv M:
 [[-1.57894737 -0.07894737  1.23684211  1.10526316]
 [-0.63157895 -0.13157895  0.39473684  0.84210526]
 [ 0.68421053  0.18421053 -0.55263158 -0.57894737]
 [ 0.52631579  0.02631579 -0.07894737 -0.36842105]]


P = np.array([[0,1,0],
              [0,0,0],
              [1,0,1]])     # P has no inverse since it is singular. 

print('det(P):', det(P))

inv(P)

det(P): 0.0

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-13-da87415bf3f7> in <module>
      5 print('det(P):', det(P))
      6 
----> 7 inv(P)

<__array_function__ internals> in inv(*args, **kwargs)

~/miniforge3/lib/python3.8/site-packages/numpy/linalg/linalg.py in inv(a)
    543     signature = 'D->D' if isComplexType(t) else 'd->d'
    544     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 545     ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    546     return wrap(ainv.astype(result_t, copy=False))
    547 

~/miniforge3/lib/python3.8/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_singular(err, flag)
     86 
     87 def _raise_linalgerror_singular(err, flag):
---> 88     raise LinAlgError("Singular matrix")
     89 
     90 def _raise_linalgerror_nonposdef(err, flag):

LinAlgError: Singular matrix


from numpy.linalg import cond, matrix_rank

A = np.array([[1,1,0],
              [0,1,0],
              [1,0,1]])

print('Condition number:\n', cond(A))

print('Rank:\n', matrix_rank(A))

y = np.array([[1, 2, 1]]).T    # need a column vector. Easier to create it using .T 

A_y = np.concatenate((A, y), axis = 1)  # Notice the axis 
print('Augmented matrix:\n', A_y)

Condition number:
 4.048917339522305
Rank:
 3
Augmented matrix:
 [[1 1 0 1]
 [0 1 0 2]
 [1 0 1 1]]


import numpy as np 
from numpy.linalg import inv

import scipy 

from scipy.linalg import lu


help(lu)

Help on function lu in module scipy.linalg.decomp_lu:

lu(a, permute_l=False, overwrite_a=False, check_finite=True)
    Compute pivoted LU decomposition of a matrix.
    
    The decomposition is::
    
        A = P L U
    
    where P is a permutation matrix, L lower triangular with unit
    diagonal elements, and U upper triangular.
    
    Parameters
    ----------
    a : (M, N) array_like
        Array to decompose
    permute_l : bool, optional
        Perform the multiplication P*L (Default: do not permute)
    overwrite_a : bool, optional
        Whether to overwrite data in a (may improve performance)
    check_finite : bool, optional
        Whether to check that the input matrix contains only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.
    
    Returns
    -------
    **(If permute_l == False)**
    
    p : (M, M) ndarray
        Permutation matrix
    l : (M, K) ndarray
        Lower triangular or trapezoidal matrix with unit diagonal.
        K = min(M, N)
    u : (K, N) ndarray
        Upper triangular or trapezoidal matrix
    
    **(If permute_l == True)**
    
    pl : (M, K) ndarray
        Permuted L matrix.
        K = min(M, N)
    u : (K, N) ndarray
        Upper triangular or trapezoidal matrix
    
    Notes
    -----
    This is a LU factorization routine written for SciPy.
    
    Examples
    --------
    >>> from scipy.linalg import lu
    >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
    >>> p, l, u = lu(A)
    >>> np.allclose(A - p @ l @ u, np.zeros((4, 4)))
    True


A = np.array([ 
              [  1,  1,   1,   1 ], 
              [  2,  3,   1,   5 ], 
              [ -1,  1,  -5,   3 ], 
              [  3,  1,   7,  -2 ]  
             ], dtype=float)


P, L, U = lu(A)

print("A:")
print(A)

print("P:")   # What is this guy? 
print(P)

print("L:")
print(L)

print("U:")
print(U)

# Is the answer the same as what we got?

A:
[[ 1.  1.  1.  1.]
 [ 2.  3.  1.  5.]
 [-1.  1. -5.  3.]
 [ 3.  1.  7. -2.]]
P:
[[0. 0. 0. 1.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [1. 0. 0. 0.]]
L:
[[ 1.          0.          0.          0.        ]
 [ 0.66666667  1.          0.          0.        ]
 [-0.33333333  0.57142857  1.          0.        ]
 [ 0.33333333  0.28571429  0.5         1.        ]]
U:
[[ 3.          1.          7.         -2.        ]
 [ 0.          2.33333333 -3.66666667  6.33333333]
 [ 0.          0.         -0.57142857 -1.28571429]
 [ 0.          0.          0.          0.5       ]]


L = P.dot(L)   # Let's just do this once


L.dot(U)     # Produces the same result though

array([[ 1.,  1.,  1.,  1.],
       [ 2.,  3.,  1.,  5.],
       [-1.,  1., -5.,  3.],
       [ 3.,  1.,  7., -2.]])

A

array([[ 1.,  1.,  1.,  1.],
       [ 2.,  3.,  1.,  5.],
       [-1.,  1., -5.,  3.],
       [ 3.,  1.,  7., -2.]])


L.dot(U) == A   # !!! the horror!

array([[ True, False, False, False],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True]])


help(np.allclose)

Help on function allclose in module numpy:

allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False)
    Returns True if two arrays are element-wise equal within a tolerance.
    
    The tolerance values are positive, typically very small numbers.  The
    relative difference (`rtol` * abs(`b`)) and the absolute difference
    `atol` are added together to compare against the absolute difference
    between `a` and `b`.
    
    NaNs are treated as equal if they are in the same place and if
    ``equal_nan=True``.  Infs are treated as equal if they are in the same
    place and of the same sign in both arrays.
    
    Parameters
    ----------
    a, b : array_like
        Input arrays to compare.
    rtol : float
        The relative tolerance parameter (see Notes).
    atol : float
        The absolute tolerance parameter (see Notes).
    equal_nan : bool
        Whether to compare NaN's as equal.  If True, NaN's in `a` will be
        considered equal to NaN's in `b` in the output array.
    
        .. versionadded:: 1.10.0
    
    Returns
    -------
    allclose : bool
        Returns True if the two arrays are equal within the given
        tolerance; False otherwise.
    
    See Also
    --------
    isclose, all, any, equal
    
    Notes
    -----
    If the following equation is element-wise True, then allclose returns
    True.
    
     absolute(`a` - `b`) <= (`atol` + `rtol` * absolute(`b`))
    
    The above equation is not symmetric in `a` and `b`, so that
    ``allclose(a, b)`` might be different from ``allclose(b, a)`` in
    some rare cases.
    
    The comparison of `a` and `b` uses standard broadcasting, which
    means that `a` and `b` need not have the same shape in order for
    ``allclose(a, b)`` to evaluate to True.  The same is true for
    `equal` but not `array_equal`.
    
    `allclose` is not defined for non-numeric data types.
    
    Examples
    --------
    >>> np.allclose([1e10,1e-7], [1.00001e10,1e-8])
    False
    >>> np.allclose([1e10,1e-8], [1.00001e10,1e-9])
    True
    >>> np.allclose([1e10,1e-8], [1.0001e10,1e-9])
    False
    >>> np.allclose([1.0, np.nan], [1.0, np.nan])
    False
    >>> np.allclose([1.0, np.nan], [1.0, np.nan], equal_nan=True)
    True


np.allclose(L.dot(U), A)  # Why?

True


y = np.array([10, 31, -2, 18])       # Is something wrong here? 
print(y)

[10 31 -2 18]


m = inv(L).dot(y)


print(m)

[18.         19.         -6.85714286  2.        ]


x = inv(U).dot(m)


print(x)   # Match this with our on-paper results.

[1. 2. 3. 4.]


yp = A.dot(x)
print(yp)

[10. 31. -2. 18.]


np.allclose(y, yp)

True


a = [
        [ 8,  3, -3 ], 
        [-2, -8,  5 ], 
        [ 3,  5,  10]
    ]

# Find diagonal coefficients
diag = np.diag(np.abs(a)) 
print(diag)

[ 8  8 10]


np.sum(np.abs(a), axis=1)

array([14, 15, 18])


# Find row sum without diagonal
off_diag = np.sum(np.abs(a), axis=1) - diag 
print(off_diag)

[6 7 8]


if np.all(diag > off_diag):
    print('matrix is diagonally dominant')    # will converge (sufficient condition)
else:
    print('NOT diagonally dominant')          # may or may not converge (because NOT a necessary condition)

matrix is diagonally dominant


# Let's apply the iterative method starting at 0 

x1 = 0
x2 = 0
x3 = 0
epsilon = 0.001   # also try better precision
converged = False

x_old = np.array([x1, x2, x3])


max_iterations = 50


print('Iteration results')
print('k      x1        x2        x3 ')

for k in range(1, max_iterations):
    
    x1 = (14 -  3*x2  +  3*x3)  /  8
    x2 = (5  +  2*x1  -  5*x3)  / (-8)
    x3 = (-8 -  3*x1  -  5*x2)  / 10
    
    x = np.array([x1, x2, x3])
    
    print("%2d   %.4f   %.4f   %.4f" % (k, x1, x2, x3))
    
    # check if it is smaller than threshold
    dx = np.sqrt(np.dot(x-x_old, x-x_old))
    
    if dx < epsilon:
        converged = True
        print('Converged!')
        break
        
    # assign the latest x value to the old value
    x_old = x

if not converged:
    print('Not converged, increase the # of iterations')

Iteration results
k      x1        x2        x3 
 1   1.7500   -1.0625   -0.7937
 2   1.8508   -1.5838   -0.5633
 3   2.1327   -1.5103   -0.6847
 4   2.0596   -1.5678   -0.6340
 5   2.1002   -1.5463   -0.6569
 6   2.0835   -1.5565   -0.6468
 7   2.0911   -1.5520   -0.6513
 8   2.0878   -1.5540   -0.6493
 9   2.0893   -1.5531   -0.6502
10   2.0886   -1.5535   -0.6498
Converged!


8*x1 + 3*x2 - 3*x3

13.997673687371334


-2*x1 - 8*x2 + 5*x3

5.001948405718945


3*x1 + 5*x2 + 10*x3

-8.0


a = [
        [ 1,  3, -3 ], 
        [-2, -1,  5 ], 
        [ 3,  5,  1]
    ]

# Find diagonal coefficients
diag = np.diag(np.abs(a)) 
print(diag)

off_diag = np.sum(np.abs(a), axis=1) - diag 
print(off_diag)

if np.all(diag > off_diag):
    print('matrix is diagonally dominant')    
else:
    print('NOT diagonally dominant')

[1 1 1]
[6 7 8]
NOT diagonally dominant


# Let's apply the iterative method starting at 0 

x1 = 0
x2 = 0
x3 = 0
epsilon = 0.001   # also try better precision
max_iterations = 10
converged = False

x_old = np.array([x1, x2, x3])

print('Iteration results')
print('k      x1        x2        x3 ')

for k in range(1, max_iterations):
    
    x1 = (14 -  3*x2  +  3*x3)  /  1
    x2 = (5  +  2*x1  -  5*x3)  / (-1)
    x3 = (-8 -  3*x1  -  5*x2)  / 1
    
    x = np.array([x1, x2, x3])
    
    print("%2d   %.4f   %.4f   %.4f" % (k, x1, x2, x3))
    
    # check if it is smaller than threshold
    dx = np.sqrt(np.dot(x-x_old, x-x_old))
    
    if dx < epsilon:
        converged = True
        print('Converged!')
        break
        
    # assign the latest x value to the old value
    x_old = x

if not converged:
    print('Not converged, increase the # of iterations')
    
# It may still have converged but it didn't ...

Iteration results
k      x1        x2        x3 
 1   14.0000   -33.0000   115.0000
 2   458.0000   -346.0000   348.0000
 3   2096.0000   -2457.0000   5989.0000
 4   25352.0000   -20764.0000   27756.0000
 5   145574.0000   -152373.0000   325135.0000
 6   1432538.0000   -1239406.0000   1899408.0000
 7   9416456.0000   -9335877.0000   18430009.0000
 8   83297672.0000   -74445304.0000   122333496.0000
 9   590336414.0000   -569005353.0000   1074017515.0000
Not converged, increase the # of iterations


# Just for reference, here's the code for crout's method


def crout(A):
    
    n = A.shape[0]
    
    U = np.zeros((n, n), dtype=np.double)
    L = np.zeros((n, n), dtype=np.double)
    
    for k in range(n):
        
        L[k, k] = A[k, k] - L[k, :] @ U[:, k]
        
        U[k, k:] = (A[k, k:] - L[k, :k] @ U[:k, k:]) / L[k, k]
        L[(k+1):, k] = (A[(k+1):, k] - L[(k+1):, :] @ U[:, k]) / U[k, k]
    
    return L, U


B = [[ 1, 2,  3],
       [ 1, 5, 20],
       [ 1,  9,  7]]
B = np.array(B)
L, U = crout(B)

L

array([[  1.        ,   0.        ,   0.        ],
       [  1.        ,   3.        ,   0.        ],
       [  1.        ,   7.        , -35.66666667]])

U

array([[1.        , 2.        , 3.        ],
       [0.        , 1.        , 5.66666667],
       [0.        , 0.        , 1.        ]])


L @ U

array([[ 1.,  2.,  3.],
       [ 1.,  5., 20.],
       [ 1.,  9.,  7.]])

Linear Algebra Basics¶

Why Study Linear Algebra?¶

Vectors¶

Vector Norm¶

Euclidean Distance¶

Manhattan Distance¶

Numpy Vectors Revisited¶

Vector Operations¶

Matrices¶

Matrix Multiplication¶

Determinants¶

Inverse, Singular, Ill-conditioned¶

Rank, Full Rank, Condition Number, Augmented Matrix¶

Solving Linear Equations¶

LU Decomposition¶

Gauss-Seidel Method¶