A system of equations

$2x+y-z=1$

$3x+4y+2z=13 \quad \quad \Rightarrow \quad \quad \left[ \begin{matrix} 2 \quad 1 -1 \\ 3 \quad 4 \quad 2 \\ 1 -5 -2 \end{matrix} \right]\left[ \begin{matrix} x \\ y \\ z \end{matrix} \right]=\left[ \begin{matrix} 1 \\ 13 \\ 0 \end{matrix} \right]$

$x-5y-2z=0 \quad \quad \quad \quad \quad \quad \quad \: \: \: \: \: \: X \quad \quad \: \beta \quad = \quad y$

$X \beta = y$

$\quad X^TX\beta = X^T y$

$\quad \beta = (X^TX)^{-1}X^T y \quad \quad $Least-squares solution

In [1]:
import numpy as np
import matplotlib.pyplot as plt
In [2]:
X = np.array([[2,1,-1],[3,4,2],[1,-5,-2]])
y = np.array([1,13,0])

beta1 = np.linalg.inv(X.T@X) @ (X.T@y)
print(beta1)

beta2 = np.linalg.solve(X,y)
print(beta2)

[ 3. -1.  4.]
[ 3. -1.  4.]

$\begin{cases} 3x_1-y_1=6 \\ -x_1+y_1=2 \end{cases}$

$\begin{cases} 3x_2-y_2=6 \\ 1.5x_2+.5y_2=3 \end{cases}$

$\begin{cases} 3x_3-y_3=6 \\ 1.5x_3-.5y_3=2 \end{cases}$

In [15]:
X1 = np.array([[3,-1],[-1,1]])
y1 = np.array([6,2])
               
X2 = np.array([[3,-1],[1.5,-.5]])
y2 = np.array([6,3])
                            
X3 = np.array([[3,-1],[1.5,-.5]])
y3 = np.array([6,2])
In [16]:
xlim = np.array([0,10])

#1)

# covert to y=mx+b format
yy1 = -X1[0,0]/X1[0,1]*xlim + y1[0]/X1[0,1] # (x = 0) y = -6, (x = 10) y = 24
print('(x = ' + str(xlim[0]) + ') y = ' + str(yy1[0]) + ', (x = ' + str(xlim[1]) + ') y = ' + str(yy1[1]))
yy2 = -X1[1,0]/X1[1,1]*xlim + y1[1]/X1[1,1] # (x = 0) y = 2, (x = 10) y = 12
print('(x = ' + str(xlim[0]) + ') y = ' + str(yy2[0]) + ', (x = ' + str(xlim[1]) + ') y = ' + str(yy2[1]))

plt.plot(xlim,yy1,label='eq 1')
plt.plot(xlim,yy2,label='eq 2')

b = np.linalg.solve(X1,y1)
print(b)
plt.plot(b[0],b[1],'o',markersize=8,label='Solution')

plt.legend()
plt.grid()
plt.show()

(x = 0) y = -6.0, (x = 10) y = 24.0
(x = 0) y = 2.0, (x = 10) y = 12.0
[4. 6.]
In [17]:
#2) Singular matrix

# covert to y=mx+b format
yy1 = -X2[0,0]/X2[0,1]*xlim + y2[0]/X2[0,1]
print('(x = ' + str(xlim[0]) + ') y = ' + str(yy1[0]) + ', (x = ' + str(xlim[1]) + ') y = ' + str(yy1[1]))
yy2 = -X2[1,0]/X2[1,1]*xlim + y2[1]/X2[1,1]
print('(x = ' + str(xlim[0]) + ') y = ' + str(yy2[0]) + ', (x = ' + str(xlim[1]) + ') y = ' + str(yy2[1]))

plt.plot(xlim,yy1,label='eq 1')
plt.plot(xlim,yy2,label='eq 2')

# b = np.linalg.solve(X2,y2) # invalid
# print(b)
# plt.plot(b[0],b[1],'o',markersize=8,label='Solution')

plt.legend()
plt.grid()
plt.show()

(x = 0) y = -6.0, (x = 10) y = 24.0
(x = 0) y = -6.0, (x = 10) y = 24.0
In [18]:
#3) Singular matrix

# covert to y=mx+b format
yy1 = -X3[0,0]/X3[0,1]*xlim + y3[0]/X3[0,1]
print('(x = ' + str(xlim[0]) + ') y = ' + str(yy1[0]) + ', (x = ' + str(xlim[1]) + ') y = ' + str(yy1[1]))
yy2 = -X3[1,0]/X3[1,1]*xlim + y3[1]/X3[1,1]
print('(x = ' + str(xlim[0]) + ') y = ' + str(yy2[0]) + ', (x = ' + str(xlim[1]) + ') y = ' + str(yy2[1]))

plt.plot(xlim,yy1,label='eq 1')
plt.plot(xlim,yy2,label='eq 2')

# b = np.linalg.solve(X2,y2) # invalid
# print(b)
# plt.plot(b[0],b[1],'o',markersize=8,label='Solution')

plt.legend()
plt.grid()
plt.show()

(x = 0) y = -6.0, (x = 10) y = 24.0
(x = 0) y = -4.0, (x = 10) y = 26.0