# Identifying special matrices

Test if a 4×4 matrix is singular, i.e. to determine if an inverse exists, before calculating it.  
Use the method of converting a matrix to echelon form, and testing if this fails by leaving zeros that can’t be removed on the leading diagonal.


```python
% matplotlib inline

import numpy as np
import numpy.linalg as la
```


```python
# Print option; suppress small numbers and exp, precision, numbers shown in matrix
# https://docs.scipy.org/doc/numpy/reference/generated/numpy.set_printoptions.html
# np.set_printoptions()  # formatter gets reset
np.set_printoptions(suppress=True, precision=3, threshold=10)
```


```python
import sympy as sp
from sympy import Matrix as spm
sp.init_printing(use_unicode=True)
from mpmath import *
mp.dps = 3 # set decimal places
```


```python
%precision 3
print(mp)
mpf(2.00000001)
2.0000001
```


```python
class MatrixIsSingular(Exception): 
    """Exception extension for errors if the matrix is singular."""
    def __init__(self, wth, msg=None):
        if msg is None:
            # Pass error message
            msg = 'Runtime error: {}'.format(wth)
        super(MatrixIsSingular, self).__init__(msg)
        self.wth = wth
```

### Manual RREF


```python
def isSingular(A):
    """Try RREF, when diagonal can't be made 1, the matrix is singular"""
    # Work with copy of mutable input
    B = np.array(A, dtype=np.float_)
    try:
        fixRowZero(B)
        fixRowOne(B)
        fixRowTwo(B)
        fixRowThree(B)
    except MatrixIsSingular:
        return True
    return False
```


```python
def fixRowZero(A):
    """"""
    if A[0,0] == 0:                 # test diagonal for non-zero
        A[0] = A[0] + A[1]          # add lower vector to populate diagonal
    if A[0,0] == 0:                 # repeat for all rows ""
        A[0] = A[0] + A[2]
    if A[0,0] == 0:
        A[0] = A[0] + A[3]
    if A[0,0] == 0:
        raise MatrixIsSingular('fixRowZero')
    A[0] = A[0] / A[0,0]            # Normalize row to set the diagonal element to one
    return A
```


```python
def fixRowOne(A):
    A[1] = A[1] - A[1,0] * A[0]     # subtract upper vector to make lower-triangle zero
    
    if A[1,1] == 0:                 # test diagonal for non-zero
        A[1] = A[1] + A[2]          # add lower vector to populate diagonal
        A[1] = A[1] - A[1,0] * A[0] # subtract upper vector to make lower-triangle zero
        
    if A[1,1] == 0:                 # again test diagonal for non-zero
        A[1] = A[1] + A[3]          # add lower vector to populate diagonal
        A[1] = A[1] - A[1,0] * A[0] # subtract upper vector to make lower-triangle zero
        
    if A[1,1] == 0:                 # final test diagonal for non-zero
        raise MatrixIsSingular('fixRowOne')
    A[1] = A[1] / A[1,1]            # Normalize row to set the diagonal element to one
    return A
```


```python
def fixRowTwo(A) :
    A[2] = A[2] - A[2,0] * A[0]
    A[2] = A[2] - A[2,1] * A[1]
    
    if A[2,2] == 0:
        A[2] = A[2] + A[3]
        A[2] = A[2] - A[2,0] * A[0]
        A[2] = A[2] - A[2,1] * A[1]
        
    if A[2,2] == 0:
        raise MatrixIsSingular('fixRowTwo')
    A[2] = A[2] / A[2,2]
    return A
```


```python
def fixRowThree(A) :
    A[3] = A[3] - A[3,0] * A[0]
    A[3] = A[3] - A[3,1] * A[1]
    A[3] = A[3] - A[3,2] * A[2]
    
    if A[3,3] == 0:
        raise MatrixIsSingular('fixRowThree')
    A[3] = A[3] / A[3,3]
    return A
```

#### Test your code


```python
A = np.array([
        [0, 7, -5, 3],
        [2, 8, 0, 4],
        [3, 12, 0, 5],
        [1, 3, 1, 3]
    ], dtype=np.float_)
isSingular(A)
spm(A)
```


```python
spm(fixRowZero(A))
## row 1 added to 0 and devided by A[0,0]
```


```python
spm(fixRowOne(A).round(2))
```


```python
spm(fixRowTwo(A).round(2))
```


```python
spm(fixRowThree(A).round(2))
```

## sympy 

http://docs.sympy.org/0.7.5/tutorial/matrices.html


```python
A = np.array([
        [2, 0, 0, 0],
        [0, 3, 0, 0],
        [0, 0, 4, 4],
        [0, 0, 5, 5]
    ], dtype=np.float_)
isSingular(A)
spm(A)
```


```python
spm(A).rref()[1]         # indices of the pivot columns
spm(A).rref()[0]         # rref
spm(A).nullspace()[0]    # nullspace
spm(A).diagonalize()     #
spm(A).det()             # determinant
if spm(A).det() == 0:
    raise MatrixIsSingular('determinant is 0')
spm(A)**-1           # matrix inverse
spm(A)**-1 @ spm(A)  # A-1 @ A = I - check inverse/rref
```

### Invertible matrix - TODO nxm matrices


```python
def is_invertible(B):
    Bp = np.linalg.inv(B)
    I = np.eye(B.shape[1], B.shape[1])
    
    assert spm(B).det != 0, 'matrix is singular'
    np.testing.assert_array_almost_equal(B @ I, B, err_msg='not invertible')
    np.testing.assert_array_almost_equal(I, Bp @ B, err_msg='not invertible')
    np.testing.assert_array_almost_equal(I, B @ Bp, err_msg='not invertible')
    np.testing.assert_array_almost_equal(Bp @ B, B @ Bp, err_msg='not invertible')
    return True

spm(B)
is_invertible(B)
```

### Transformation of base


```python
A = np.array([
        [2, 1],
        [1, 2]
    ], dtype=np.float_)
```


```python
# New base
B = np.array([
        [3, 1],
        [1, 1]
    ], dtype=np.float_)
B
```


```python
# Rotation 45
R = np.array([
        [1, -1],
        [1, 1]
    ], dtype=np.float_) #/np.sqrt(2)
R
```


```python
# Identity
I = np.eye(2)
I
```


```python
# Back to base
Binv = np.linalg.inv(B)
```


```python
spm(A)
spm((Binv @ I @ B).round(2))
spm(Binv @ R @ B)
spm((Binv @ A @ B).round(2))
```

## Eigen experiments

Eigen vectors `x` and eigen values `lambda`:
```python
A @ x = lambda @ x
A @ x - lambda @ I @ x = 0
```


```python
%precision 3
B_eig = spm(B).eigenvects()[0][2][0]
B_eigv = spm(B).eigenvects()[0][0]
spm(B).eigenvects()
B @ B_eig
B_eigv * B_eig
(B @ B_eig - B_eigv @ spm(np.eye(B.shape[0])) @ B_eig)

l, v = la.eig(B)
l
v
B
l[0] * v[:,0], B @ v[:,0]
np.testing.assert_array_almost_equal(l[0] * v[:,0], B @ v[:,0])


```

### SVD Singular Value Decomposition
When a is a 2D array, it is factorized as u @ np.diag(s) @ vh = (u * s) @ vh, 
where u and vh are 2D unitary arrays and s is a 1D array of a‘s singular values. 


```python
u, s, vh = np.linalg.svd(B, full_matrices=True)
u
s
vh
u @ s @ vh
```


```python
# Transform in original
a, b, c, d = 1, 0, 2, -1
T = np.array([
        [a, b],
        [c, d]
    ], dtype=np.float_)
```


```python
# e = values, v = unit eigen matrix (normalized)
e, C = la.eig(T)
'C'
spm(C.round(2))
'C_inv'
spm(la.inv(C).round(2))
'Eigenvalues: {}'.format(e)

```


```python
%precision 3
# Change of basis matrix with Eigen vectors of T
a, b, c, d = 1, 2, 0, 1
C = np.array([
        [a, b],
        [c, d]
    ], dtype=np.float_)
'C'
spm(C)

# Inverse Eigen matrix
Cinv = la.inv(C)
'Cinv'
spm(Cinv)

# Diagonalized matrix of eigen values 
D = la.inv(C) @ T @ C
'D'
spm(D)

A = C @ D @ la.inv(C)
'A'
spm(A)

An = C @ D @ D @ D @ D @ D @ la.inv(C)
'An'
spm(An)

```


```python
# Transform in original
a, b, c, d, e, f, g, h, i = 4,-5,6, 7,-8,6, 3/2,-1/2,-2
A = np.array([
        [a, b, c],
        [d, e, f],
        [g, h, i]
    ], dtype=np.float_)
vals, vecs = la.eig(A)
print(vals)
print(vecs)
print('sqrt:', 6, 1/np.sqrt(6))

# Normalized
i, j, k = 0.5,-.5,-1
v = np.array([[i, j, k]])
print('normalized: ', v/la.norm(v))
```

### Eigen vector - experiment
```python
Eigen vectors A@x = lambda@x => (A@x - lambda@I@x) = 0
e = values, v = unit eigen matrix (normalized)
```

#### Complex to real numbers


```python
A = np.array([[0,0,0,1], [1,0,0,0], [0,1,0,0], [0,0,1,0]])
vals, vecs = la.eig(A)
vals
# np.real() converts complex numbers to real numbers
np.real(vals)
np.real(vecs)
```


```python
A = np.array([[0.0,0.0,0.0,1.0], [1.0,0.0,0.0,0.0], [0.0,1.0,0.0,0.0], [0.0,0.0,1.0,0.0]])
D = np.array([[0.1,0.1,0.1,0.7], [0.7,0.1,0.1,0.1], [0.1,0.7,0.1,0.1], [0.1,0.1,0.7,0.1]])
d = 0.5

M = d * A + (d - 1) * D/4
spm(A), spm(D), spm(M.round(2))
vals, vecs = la.eig(M)

np.real(vals)
np.real(vecs)
```


```python
d = 0.9
M = d * A + (d - 1) * D/4
vals, vecs = la.eig(M)

np.real(vals)
spm(np.real(vecs).round(2))
```


```python
a, b, c, d = 3/2, -1, -1/2, 1/2
T = np.array([
        [a, b],
        [c, d]
    ], dtype=np.float_)
spm(T)
```


```python

vals, vecs = la.eig(T)
print('vals', vals)
C = vecs
spm(C)

D = la.inv(C) @ T @ C
spm(D)

D_2 = vals * np.eye(2)
spm(D_2)

A2 = C @ D @ D @ la.inv(C)
spm(A2)

spm(T @ T)
```

### Eigen stuff - sympy vs. numpy.linalg


```python
-1 + sp.sqrt(3)
-1 + np.sqrt(3)
```


```python
a, b, c, d, e, f, g, h, i = 4,-5,6, 7,-8,6, 3/2,-1/2,-2
A = np.array([
        [a, b, c],
        [d, e, f],
        [g, h, i]
    ], dtype=np.float_)
spm(A)
```


```python
spm(np.array([4,-5,6, 7,-8,6, 3/2,-1/2,-2], dtype=np.float_).reshape(-1,3))
```


```python
# sympy - algebraic('non-normalised') eigenvectors
# numpy.linalg - numeric eigenvectors
sorted(spm(A).eigenvects(), key=lambda x: -x[0])
```

### Eigen functions


```python
def eigen(X, descending=True):
    """Eigen values and -vectors in descending order"""
    from collections import defaultdict
    eigen_dict = defaultdict()
    for i, v in enumerate(spm(A).eigenvects()):
        eigen_dict[spm(X).eigenvects()[i][0]] = spm(X).eigenvects()[i][2][0]
    return sorted(eigen_dict.items(), key=lambda x: [1, -1][descending] * x[0])

def eigen_sym(X, normalised=False, descending=True):
    """Eigen values and -vectors in descending order"""
    from collections import defaultdict
    eigen_dict = defaultdict()
    for i, v in enumerate(spm(X).eigenvects()):
        norm = [1, (spm(X).eigenvects()[i][2][0].T @ spm(X).eigenvects()[i][2][0])**.5][normalised]
        eigen_dict[spm(X).eigenvects()[i][0]] = spm(X).eigenvects()[i][2][0]/norm
    return sorted(eigen_dict.items(), key=lambda x: [1, -1][descending] * x[0])

def eigen_la(X, descending=True):
    """Eigen values and -vectors in descending order"""
    vals, vecs = la.eig(X)
    idx = np.argsort([1, -1][descending] * vals)
    return [(val.round(), spm(vec)) for val, vec in zip(vals[idx], vecs.T[idx])]
```


```python
eigen_sym(A)
eigen_sym(A, False)
eigen_sym(A, True)
```


```python
eigen_la(A, True)
```

### Angle (in radius) between vectors x and z


```python
def v_len(x):
    return (x.T @ x)**.5

def radius(x, y):
    x_, y_ = v_len(x), v_len(y)
    if x_ == y_: return 0
    return np.arccos(x.T @ y / (x_ * y_))
```


```python
x = np.array([1,0,6])
y = np.array([-1,0,8])
z = x - y

radius(x, y)
```


```python

```

2da6e34c50dbced44dbfc78f2587243c31765324
ipynb
Jupyter Notebook
_math/Imperial College London - Math for ML/matrix_rref_eigen_sympy_errorclass.ipynb
aixpact/data-science
f04a54595fbc2d797918d450b979fd4c2eabac15
MIT
2020-07-22T23:12:39.000Z
2020-07-25T02:30:48.000Z
_math/Imperial College London - Math for ML/matrix_rref_eigen_sympy_errorclass.ipynb
aixpact/data-science
f04a54595fbc2d797918d450b979fd4c2eabac15
MIT
_math/Imperial College London - Math for ML/matrix_rref_eigen_sympy_errorclass.ipynb
aixpact/data-science
f04a54595fbc2d797918d450b979fd4c2eabac15
MIT
Qwen/Qwen-72B
1. YES
2. YES