Linear Algebra with Python Libraries - Numpy and matplotlib
Bikash Santra and Avisek Gupta
Indian Statistical Institute, Kolkata
1. Linear Independence of vectors
# space spanned by v1 and v2 shown in yellow
%matplotlib inline
import numpy as np
v1 = np.array([1,2])
v2 = np.array([2,4])
import matplotlib as mpl
import matplotlib.pyplot as plt
fig = plt.figure()
maxlim = max(v1.max(), v2.max())
minlim = max(v1.min(), v2.min())
lim = max(np.abs(maxlim), np.abs(minlim))
points = np.linspace(-lim, lim, 50)
p0 = np.repeat(points, 50)
p1 = np.array(list(points) * 50)
space = np.hstack((p0[:,None], p1[:,None]))
mat = np.hstack((v1[:,None],v2[:,None]))
mask = np.zeros((space.shape[0]), dtype=np.int)
for apoint in range(space.shape[0]):
if np.linalg.matrix_rank( np.hstack((mat, space[apoint,:][:,None])) ) == np.linalg.matrix_rank(mat):
mask[apoint] = 1
plt.scatter(space[mask==1,0], space[mask==1,1], c='#FFFF00',s=60)
plt.plot([0, 0], [-lim, lim], c='#dddddd')
plt.plot([-lim, lim], [0, 0], c='#dddddd')
plt.plot([0, -v2[0]], [0, -v2[1]], linewidth=4, c='#FFFF00')
plt.plot([0, -v1[0]], [0, -v1[1]], linewidth=4, c='#FFFF00')
plt.plot([0, v2[0]], [0, v2[1]], linewidth=4, c='r')
plt.plot([0, v1[0]], [0, v1[1]], linewidth=4, c='b')
plt.xlim(-lim, lim)
plt.ylim(-lim, lim)
plt.xlabel('x', size=16)
plt.ylabel('y', size=16)
plt.show()
2. Matrix ranks
A = np.reshape(np.arange(1,10), (3,3)).T
print('A =\n', A)
print('rank(A) =', np.linalg.matrix_rank(A))
print('')
B = np.array([[1, 4, 2], [2, 5, 1], [3, 7, 4]])
print('B = \n', B)
print('rank(B) =', np.linalg.matrix_rank(B))
3. Eigenvalues and eigenvectors
B = np.array([[1, 4, 2], [2, 5, 1], [3, 7, 4]])
print('B = \n', B)
eigvals, eigvecs = np.linalg.eigh(B)
print('eigenvalues of B =', eigvals)
print('eigenvectors of B =\n', eigvecs)
4. Solving Linear Equations
B = np.array([[1, 4, 2], [2, 5, 1], [3, 7, 4]])
b = np.array([2,4,6])
print(np.linalg.solve(B, b))
We can use $x = A^{-1}b$ to solve a system of linear equation $Ax=b$ only when $A$ is a square matrix with linearly independent columns (and rows).
We cannot use $x = A^{-1}b$ solve $Ax=b$ when :
- $A$ is square, but does not have linearly independent columns (and rows).
- $A$ is not square.
A = np.reshape(np.arange(1,10), (3,3)).T
b = np.array([1,2,3])
# Will produce an error since A is singular
#print(np.linalg.solve(A, b))
# Cannot find inverse of non-square matrices
#C = np.random.rand(3,4)
#print(np.linalg.inv(C))
We can always find a solution (if it exists) to $Ax=b$ using $A^{\dagger}$, which is the pseudo-inverse of $A$, from the following equation:
$x = A^{\dagger}b$
Here $b$ is projected to the column space of $A$, to find a solution $x$.
$A$ need not have linearly-independent columns, and it can be rectangular. $b$ will be projected to the span of the columns of $A$ to find $x$.
# pseudo-inverse to solve Ax=b
A = np.reshape(np.arange(1,10), (3,3)).T
print(A)
b = np.array([1,3,5])
pA = np.linalg.pinv(A)
x = np.dot(pA, b)
print(x)
# pseudo-inverse to solve Ax=b when A is rectangular
A = np.array([[1, 4, 2], [2, 5, 1], [3, 7, 4], [5, 2, 3]])
print(A)
b = np.array([2,4,6,3])
pA = np.linalg.pinv(A)
x = np.dot(pA, b)
print(x)
# solving Ax=b for random A and b
n = np.random.randint(3,5+1)
d = np.random.randint(3,5+1)
A = np.random.rand(n,d)
b = np.random.rand(n)
print('Dimensions of A =(', n,',',d,')')
print('A =', A)
print('b =', b)
pA = np.linalg.pinv(A)
x = np.dot(pA, b)
print('x =', x)
5. Linear Regression
We wish to fit a line $w^T x + b$ through some given data $(x_1,y_1), (x_2,y_2), ... ,(x_n,y_n)$.
data = np.array([
[0.20, 0.30], [0.32, 0.23], [0.39, 0.35], [0.46, 0.42],
[0.51, 0.39], [0.60, 0.47], [0.64, 0.51], [0.75, 0.46]
])
plt.figure()
plt.scatter(data[:,0], data[:,1], marker='o', c='gray')
plt.show()
The equation of a line is $y=mx+c$. To 'best' fit line through the given points might be one for which the distance to all points $(x_i,y_i)$ is minimized.
An error (or cost) function can be designed $J = \sum_{i=1}^{n} ((m x_i + c) - y_i)^2 $.
The values of $m$ and $b$ that minimize $J$ can be computed by equating the corresponding derivatives of $J$ to zero:
$\frac{\partial}{\partial m} J = 0$, and
$\frac{\partial}{\partial b} J = 0$.
Solving for $m$ and $b$:
$m = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{n}(x_i - \bar{x})^2}$
$c = \bar{y} - m\bar{x}$
# Linear Regression
%matplotlib inline
data = np.array([
[0.20, 0.30], [0.32, 0.23], [0.39, 0.35], [0.46, 0.42],
[0.51, 0.39], [0.60, 0.47], [0.64, 0.51], [0.75, 0.46]
])
# Solving the optimization criterion to get m and c
m = (data[:,0] - data[:,0].mean()).dot(data[:,1] - data[:,1].mean()) \
/ ((data[:,0] - data[:,0].mean()) ** 2).sum()
c = data[:,1].mean() - m * data[:,0].mean()
print('c =', c,', m =', m, )
plt.figure()
plt.scatter(data[:,0], data[:,1], marker='o', c='gray')
plt.plot(
np.array([data[:,0].min()-0.2,data[:,0].max()+0.2]),
m * np.array([data[:,0].min()-0.2,data[:,0].max()+0.2]) + c, c='r'
)
plt.show()
6. Linear Regression:
Rewriting the objective as: $\min ||X \theta - y||^2$
where the first column of $X$ is all $1$s.
# Linear Regression
%matplotlib inline
data = np.array([
[0.20, 0.30], [0.32, 0.23], [0.39, 0.35], [0.46, 0.42],
[0.51, 0.39], [0.60, 0.47], [0.64, 0.51], [0.75, 0.46]
])
X = np.hstack((np.ones((data.shape[0], 1)), np.reshape(data[:,0], (data.shape[0], 1))))
y = data[:,1]
theta, residual_sum, rank, s = np.linalg.lstsq(X, y, rcond=1e-9)
print(theta)
plt.figure()
plt.scatter(data[:,0], data[:,1], marker='o', c='gray')
plt.plot(
np.array([data[:,0].min()-0.2,data[:,0].max()+0.2]),
theta[1] * np.array([data[:,0].min()-0.2,data[:,0].max()+0.2]) + theta[0], c='r'
)
plt.show()