Linear Algebra in Python 1 — 方程组的几何解释

Linear Algebra in Python 1

下面是四个知识点，其中Column Picture尤其重要。

• n linear equations, n unknowns.
• Row picture
• Column picture
• Matrix form

1. Row Picture（行向量视角）

对于下面方程组，两个方程两个未知数。

\begin{align*} 2x - y &= 0 \\ -x + 2y &= 3 \end{align*}

如果我们把这两个方程看成是平面直角坐标系的两条直线的话（看成Row Picture），那么它们的图形如下（这里使用IPython Notebook）：

%matplotlib
import matplotlib.pyplot as plt
import numpy as np

x1 = np.linspace(-3, 3, 50)
y1 = 2 * x1
x2 = x1
y2 = (3 + x2)/2
fig, ax = plt.subplots()
plt.plot(x1, y1,'-', color='C1')
plt.plot(x2, y2, '-', color='C1')
plt.text(1, 2, '(1, 2)', color='C2')
plt.show()

这里可以看出两条直线的交点为$(1, 2)$点，这里正是满足方程的解：$x=1, y=2$。

2. Column Picture（列向量视角）

如果我们用列向量（Column）形式表示方程组的话可以写成下面形式：

\begin{align*} x \begin{bmatrix} 2 \\ -1 \end{bmatrix} + y \begin{bmatrix} -1 \\ 2 \end{bmatrix} {=} \begin{bmatrix} 0 \\ 3 \end{bmatrix} \end{align*}

这里有一句经典的话：

Linear Algebra is the linear combination of Columns.

线性代数就是列向量的各种（线性）组合。

接下来我们在坐标系中画出两个向量(分别记为$\vec{v1} , \vec{v2} $)。

vecs = np.array([[0, 0, -1, 2], [0, 0, 2, -1]])
X, Y, U, V = zip(*vecs)

plt.figure()
ax = plt.gca()
ax.quiver(X, Y, U, V, angles='xy', scale_units='xy', scale=1)

ax.text(-1, 2, "v1", style='italic',
            bbox={'facecolor':'red', 'alpha':0.5, 'pad':0.5})
ax.text(2, -1, "v2", style='italic',
            bbox={'facecolor':'red', 'alpha':0.5, 'pad':0.5})
            
ax.set_xticks(np.arange(vecs.min(), vecs.max() + 1, 1.0))
ax.set_yticks(np.arange(vecs.min(), vecs.max() + 1, 1.0))
ax.set_xlim([vecs.min(), vecs.max() + 1])
ax.set_ylim([vecs.min(), vecs.max() + 1])

plt.grid(b=True, which='major') #<-- plot grid lines
plt.draw()
plt.show()

如果我们对这两个向量进行向量间的操作（向量加法，因为我们上面已经知道了$x=1, y= 2$是方程的解），会得到下面的结果：

\begin{align*} 1 \begin{bmatrix} 2 \\ -1 \end{bmatrix} + 2 \begin{bmatrix} -1 \\ 2 \end{bmatrix} {=} \begin{bmatrix} 0 \\ 3 \end{bmatrix} \end{align*}

在坐标系中表示的话即是：$\vec{v1} + 2\vec{v2}$

vecs = np.array([[0, 0, -1, 2], [0, 0, 2, -1], [2, -1, -2, 4]])
X, Y, U, V = zip(*vecs)

plt.figure()
ax = plt.gca()
ax.quiver(X, Y, U, V, angles='xy', scale_units='xy', scale=1)

ax.text(-1, 2, "v1", style='italic',
            bbox={'facecolor':'red', 'alpha':0.5, 'pad':0.5})
ax.text(2, -1, "v2", style='italic',
            bbox={'facecolor':'red', 'alpha':0.5, 'pad':0.5})

ax.text(0, 3, "2v1 + v2", style='italic',
            bbox={'facecolor':'red', 'alpha':0.5, 'pad':0.5})
            
ax.set_xticks(np.arange(vecs.min(), vecs.max() + 1, 1.0))
ax.set_yticks(np.arange(vecs.min(), vecs.max() + 1, 1.0))
ax.set_xlim([vecs.min(), vecs.max() + 1])
ax.set_ylim([vecs.min(), vecs.max() + 1])

plt.grid(b=True, which='major') #<-- plot grid lines
plt.draw()
plt.show()

当然，上面我们是在知道了方程的解的情况下用$\vec{v1} + 2\vec{v2}$，然后就可以得到结果。

如果我们还不知道解的话，我们的题目就变成了寻找各种$x$, $y$组合，使

\begin{align*} x \vec{v1} + y\vec{v2}= \begin{bmatrix} 0 \\ 3 \end{bmatrix} \end{align*}

成立，我们可以想象的是，如果$\vec{v1}, \vec{v2}$可以表示平面内的所有向量（相互独立），那么它们一定可以表示$\begin{bmatrix}0 \\ 3 \end{bmatrix}$ 这个向量。

3. Matrix form（矩阵表示）

上面的方程组：

\begin{align*} 2x - y &= 0 \\ -x + 2y &= 3 \end{align*}

如果使用矩阵（Matrix）形式来表示，可以写成我们习惯的\matrix{A} \bf{x} {=} \bf{b}形式（其中A是系数矩阵，x是未知数向量，b也是一个向量）。

\begin{align*} \begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} {=} \begin{bmatrix} 0 \\ 3 \end{bmatrix} \end{align*}

从上面我们已经知道该方程组有解，称矩阵\matrix{A} 是一个非奇异、可逆矩阵（no-singular, invertible matrix），当然，如果方程不可解，称\matrix{A} 为一个奇异、不可逆矩阵（singular, not invertible matrix）。

下面我们考虑3个方程3个未知数的情况。

\begin{align*} 2x - y &= 0 \\ -x + 2y - z&= -1 \\ -3y + 4z &= 4 \end{align*}

首先，我们使用Row Picture视角，在三维坐标系中画出三个平面，它们交点即是方程组的解。

from mpl_toolkits.mplot3d import axes3d
# create x,y
xx, yy = np.meshgrid(range(10), range(10))

# calculate corresponding z
z1 = xx * 0
z2 = (-1 * xx + 2 * yy + 1) * 1.0
z3 = (3 * yy + 2 + 4) * 1.0 / 4.0

# plot the surface
plt3d = plt.figure().gca(projection='3d')
plt3d.plot_surface(xx, yy, z1, color='C1')
plt3d.plot_surface(xx, yy, z2, color='C2')
plt3d.plot_surface(xx, yy, z3, color='C3')
plt.show()

不过，老实说，交点很难看出来它们的交点（知道答案在来看就会容易点，毕竟马后炮 😃 ）。将方程写为矩阵形式：

\begin{align*} \begin{bmatrix} 2 & -1 &0 \\ -1 & 2 & -1 \\ 0 & -3 & 4 \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix} {=} \begin{bmatrix} 0 \\ -1 \\ 4 \end{bmatrix} \end{align*}

即：

\begin{align*} x \begin{bmatrix} 2 \\ -1 \\ 0 \end{bmatrix} + y \begin{bmatrix} 2 \\ -1 \\ 0 \end{bmatrix} + z \begin{bmatrix} 0 \\ -1 \\ 4 \end{bmatrix} {=} \begin{bmatrix} 0 \\ -1 \\ 4 \end{bmatrix} \end{align*}

从上面的式子中我们可以清楚地看到当$x=0, y=0, z=1$时式子成立（因为这个方程组是故意设计的），现在我们切换到Column Pitcure视角。

from mpl_toolkits.mplot3d import axes3d
vecs = np.array([[0, 0, 0, 2, -2, 0], [0, 0, 0, -1, 2, -3], [0, 0, 0, 0, -1, 4]])
X, Y, Z, U, V, W = zip(*vecs)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_xlim3d(-5, 5)
ax.set_ylim3d(-5, 5)
ax.set_zlim3d(-5, 5)
ax.quiver(X, Y, Z, U, V, W, length = 5, normalize = True)
plt.show()

这样我们可以更清楚的看到$x=0, y=0, z=1$时满足条件。

上面我们举了的例子是有解的情况，那么

Can I save Ax=b for every b? （对于每一个b，方程是否都有解？）

or do the linear combinations of the columns fill 3-D space ？（换句话说，A中的Columns通过线性组合是否能够填充满整个三维空间？

下一节将通过elimination（消元法）来解决这个问题。

4. Python求解方程

这里使用的是sympy包

• 方程一

  \begin{align*} 2x - y &= 0 \\ -x + 2y &= 3 \end{align*}

  from sympy import solve, symbols
  x, y = symbols('x y')
  sol = solve([2*x - y, -x + 2*y - 3], dict=True)
  print(sol)
  # ouput: [{y: 2, x: 1}]

• 方程二

\begin{align*} 2x - y &= 0 \\ -x + 2y - z&= -1 \\ -3y + 4z &= 4 \end{align*}

from sympy import solve, symbols
x, y, z = symbols('x y z')
sol = solve([2*x - y, -x + 2*y - z + 1, -3*y + 4*z - 4], dict=True)
print(sol)
# output: [{y: 0, x: 0, z: 1}]

5. 参考文献

方程组的几何解释

mplot3d

Solvers