本系列用彩色的 Latex 公式呈现总结 MIT Gilbert Strang 教授经典的线性代数课程18.06之精华，便于大家回顾复习。

Gauss-Jordan 算法是一种手工计算矩阵逆的技术。现在，你几乎不应该手动计算矩阵逆，即使是在计算机上，但 Gauss-Jordan 仍然有用，因为

• 它有助于我们理解逆矩阵何时以及为何存在。

• 它给了我们另一个例子来帮助我们理解消元的结构

回顾：矩阵的逆

线性算子的 \(A^{-1}\) 是“撤消”操作的操作 \(A\) ，对于任何 \(x\)，有:

即 \[\boxed{Ax=b \implies x = A^{-1} b}\]

这意味着

• 对于m×m 方阵 \(A\)， \(A^{-1}\) 存在仅当 \(A\) 具有 m 个主元的

因为对于非方阵或具有一个或多个“零主元”的矩阵，我们不能总是解决 \(Ax=b\) （会在反向替换期间除以零）。也很容易看出 \(\boxed{(A^{-1})^{-1} = A}\) , 即 \(A\) 撤消操作 \(A^{-1}\) .

等价地

这里， \(I\) 是 m×m 单位矩阵——在线性代数中，我们通常从上下文来推断 \(I\) 的形状，但如果它是模棱两可时，可以写成 \(I_m\).

乘积的逆： (AB)⁻¹ = B⁻¹A⁻¹

很容易看出乘积 \(AB\) 的逆是两者逆的反向乘积：

\[\boxed{(AB)^{-1} = B^{-1} A^{-1}}\]

直观理解为，当反转一系列操作时，总是需要以倒序的顺序回溯操作。

容易证明这个结论，因为

例如，我们看到高斯消元对应于因式分解 \(A = LU\)， 此时，\(U\) 是消除的结果，并且 \(L\) 只是消除步骤的记录。

然后

注意：逆允许我们“除以矩阵”，但我们总是要清楚我们是在左边还是右边。以下符号很方便，可用于 Julia 和 Matlab 中

通过线性方程求逆

方程 \(A A^{-1} = I\) 实际上给了我们计算的算法 \(A^{-1}\) .

假设我们将 \(A\) 以及 \(I\) 表示成列的形式，即

那么

这里的关键事实是，将 \(A\) 乘以右侧的矩阵相当于将 \(A\) 乘以该矩阵的每一列，通过写出计算可以很容易地看到这一点。

如此，有了 m 个方程组，第 k 个方程组对应 \(A x_k = e_k\)，解得这个方程组就可以得到 \(A^{-1}\) 的第 k 列 。也就是说，要找到 m×m 矩阵 \(A\) 的 \(A^{-1}\)，我们必须求解 m 个 \(Ax=b\)。

• 换句话说，对于任何矩阵 \(B\) , \(Be_k\) 为 \(B\) 的第 k 列。 所以 \(A^{-1}\) 第 k 列是 \(x_k = A^{-1} e_k\) ，即 \(Ax_k = e_k\) 的解

• 理想情况下，当我们进行一次高斯消元 \(A=LU\) 后，就可以计算出 \(x_k = U^{-1} L^{-1} e_k\) ，做法是通过在 \(I\) 的每一列上做向前和向后替换（这本质上是计算机的做法）。

举例计算 L⁻¹ = E

例如，我们如何计算我们在上一课中从高斯消元中得到的 \(L\) 矩阵的逆矩阵，我们知道，结果应该是 \(L^{-1} = E\)

对于每个 \(e_1,e_2,e_3\)，我们解对应方程组

\(e_k\) 是 3×3 \(I\) 的第 k 列。

例如对于 \(e_1\), 找到第一列 \(x_1\) 的 \(L^{-1} = E\):

The Gauss–Jordan 算法

Gauss-Jordan 可以被视为求解的一种技巧（主要用于手工计算）\(A b_k = e_k\). 但是用代数的方式思考也很好——这是我们高斯消元的“矩阵观点”的一个很好的应用。

简而言之，Gauss-Jordan 的想法是：如果我们对 \(A\) 执行一些行操作以获得 \(I\)，那么对 \(I\) 执行相同的行操作会得到 \(A^{-1}\)。为什么？

• 行操作对应于从 \(A\) 左边乘以一组矩阵 \(E=\cdots E_2 E_1\)

• 所以，对 \(A\) 做行操作将其变成 \(I\) 意思等价于 \(EA = I\)， 因此 \(E = A^{-1}\) .

• 对 \(I\) 执行相同的行操作，相当于 \(I\) 左乘矩阵 \(E\) , 即 \(EI\)，因为 \(EI = E\) 并且 \(E = A^{-1}\)，所以结果就是 \(A^{-1}\)

这就是我们可以用扩展矩阵来进行高斯消除，对 \(A\) 和 \(I\) 同时执行相同的行操作，即

\(A \to I\) 消除步骤

下面是具体将 \(A\) 通过行操作变换成 \(I\) 的步骤

1. 做普通高斯消元“向下”将 \(A\) 转变成 \(U\) （上三角矩阵）。
2. 然后，做高斯消去“向上” \(U\) 消除对角线以上的元素，将 \(U\) 转换成对角矩阵 \(D\)
3. 最后，将 \(D\) 的每一行除以对角线元素值，得到 \(I\)

Gauss-Jordan 示例

让我们执行这些 \(A \to I\) 消除步骤 \(3 \times 3\) 矩阵 \(A\) : 先向下消除获得 \(U\) , 然后向上消除得到 \(D\) ，最后除以对角线元素值得到 \(I\):

没问题！很容易看出，只要 \(A\) 具有所有主元（即它是非奇异的），这将起作用。

永远不要计算矩阵逆！

矩阵逆很有趣，它在分析操作中非常方便，因为它们允许您轻松地将矩阵从方程的一侧移动到另一侧。但是，在“严肃的”数值计算中几乎从不计算逆矩阵。每当你看到 \(A^{-1} B\) 或者 \(A^{-1} b\)，当您在计算机上实现它时，您应该阅读 \(A^{-1} B\) 作为“解决 \(AX = B\) 通过某种方法。”例如，通过 \(A \backslash B\) 或通过首先计算 \(LU\) 分解来解决它 \(A\) 然后用它来解决 \(AX = B\) .

你通常不计算逆矩阵的一个原因是它很浪费：一旦你有 \(A=LU\)（稍后我们将把它概括为“\(PA = LU\)")，你可以解决 \(AX=B\) 直接不用找 \(A^{-1}\), 和计算 \(A^{-1}\) 如果您只需要解决几个右手边，则需要更多的工作。

另一个原因是，对于很多特殊的矩阵，有办法解决 \(AX=B\) 比你能找到的要快得多 \(A^{-1}\). 例如，实践中的许多大矩阵是稀疏的（大部分为零），通常对于稀疏矩阵，您可以安排 \(L\) 和 \(U\) 也要稀疏。稀疏矩阵比一般的“密集”矩阵更有效，因为您不必乘以（甚至存储）零。即使 \(A\) 是稀疏的，然而， \(A^{-1}\) 通常是非稀疏的，因此如果计算逆矩阵，则会失去稀疏的特殊效率。

例如：

• 如果你看到 \(U^{-1} b\) 在哪里 \(U\) 是上三角，不用计算 \(U^{-1}\) 明确！只要解决 \(Ux = b\) 通过反向替换（从底行向上）。

• 如果你看到 \(L^{-1} b\) 在哪里 \(L\) 是下三角，不用计算 \(L^{-1}\) 明确！只要解决 \(Lx = b\) 通过前向替换（从顶行向下）。

本系列用彩色的 Latex 公式呈现总结 MIT Gilbert Strang 教授经典的线性代数课程18.06之精华，便于大家回顾复习。

本文在矩阵乘法和解方程组的基础上引入消元对应的矩阵乘法和矩阵的逆。

消元操作的矩阵表示

回顾上一节，解方程组的意义和过程 中，在消元步骤时，将矩阵 A 依行顺序向下，每一行通过减去上面行的若干倍，将矩阵 A 逐行规范成上三角矩阵形式。

如此，第一个矩阵 \([1 \quad 0 \quad 0]\) 乘以 \(A\) 形成了乘积矩阵的第一行，具体来说， \([1 \quad 0 \quad 0]\) 表示取 \(A\) 的1个第一行，0个第二行，0个第三行的和组成乘积矩阵第一行。 同样地，乘积矩阵第二行由第二个矩阵 \([-3 \quad 1 \quad 0]\) 形成，表示取 \(A\) 的-3个第一行，1个第二行，0个第三行之和。

具体转换如下

至此，操作 \(E_{21}\) 可以用矩阵 \(E_{21}\) 右乘以 \(A\) 来表示，即

上述矩阵乘法可以用简单地 Python 代码来核实

1
2
3
4
5
6
7
8
9
10
11
12
13
14

import numpy as np

A = np.array([
        [1., 2., 1.], 
        [3., 8., 1.],
        [0., 4., 1.]
])
e_21 = np.array([
    [1., 0., 0.], 
    [-3., 1., 0.],
    [0.,0.,1.]
])

e_21 @ A

运行结果为

1
2
3

array([[ 1.,  2.,  1.],
       [ 0.,  2., -2.],
       [ 0.,  4.,  1.]])

将 \(A\) 矩阵消元的第二个操作是 \(E_{32}\)，第三行减去2倍第二行，这个操作依样画葫芦可以写成

即 \(A\) 在两次左乘相应矩阵后变成上三角矩阵 \(U\)

由于乘法结合律，可以将两次操作 \(E_{32}, E_{21}\) 视为统一的操作 \(E\)

注意上式中， \(E_{32},E_{21}\) 为下三角矩阵。

我们再用代码来验证下矩阵乘法结合律

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

import numpy as np

e_32 = np.array([
    [1., 0., 0.], 
    [0., 1., 0.],
    [0., -2., 1.]
])

e_21 = np.array([
    [1., 0., 0.], 
    [-3., 1., 0.],
    [0.,0.,1.]
])

A = np.array([
        [1., 2., 1.], 
        [3., 8., 1.],
        [0., 4., 1.]
])

1
2

e_32 @ e_21 @ A
e_32 @ (e_21 @ A)

结果都为

1
2
3

array([[ 1.,  2.,  1.],
       [ 0.,  2., -2.],
       [ 0.,  0.,  5.]])

矩阵的逆

以基本行操作 \(E_{21}^{-1}\) 为例，此变换是从第二行中减去三倍的第一行

那么其逆变换就是给第二行加上三倍的第一行，所以逆矩阵就是

\(E_{32}^{-1}\) 也是同样

我们已经知道这两个行变换的综合效果为 \(E\)

发现综合结果的逆以相反的顺序相乘，这个也容易理解，因为逆操作将一个上三角矩阵 \(U\) 恢复成 \(A\)，这个过程是自下而上的行变换。

最后，来核实 \(E E^{-1} = I\)。

本文继续以彩色Latex 方式总结了方程组的行视角，列视角的几何意义；并回顾了解方程的两个步骤：消元和回代。内容对应于MIT 18.06 Gilbert Strang 线性代数视频课程第一，二节。

本系列链接如下

01 - 矩阵乘法的五种理解

02 - 解方程组的意义和过程

方程组两种几何解释

二元方程组

来看一个具体的二元线性方程组

\[ \begin{cases} 2x&-y&=0 \\ -x&+2y&=3 \end{cases} \]

写成矩阵形式

\[ \begin{bmatrix} 2&-1 \\ -1&2 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 0 \\ 3 \end{bmatrix} \]

行视角

回顾 $2x - y = 0 $ 为所有满足此条件的 (x, y) 的集合，即集合组成二维平面的一条直线，如下图蓝线所示。

$ -x + 2y = 3$ 则对应绿线。

因此方程组的解 x = 1, y = 2 为两条直线的交点，这就是方程组的行视角：将系数矩阵按行切分，则每一行表示一个约束条件，几何意义是N维空间的一个子空间。在二元方程中，一行表示一条线，三元方程中，一行表示一个平面（详见后一小节）。 \[ \begin{bmatrix} \color{blue} 2 x & {\color{blue} -1 y} \\ \color{green} -1 x & \color{green} 2 y \end{bmatrix} = \begin{bmatrix} \color{blue} 0 \\ \color{green} 3 \end{bmatrix} \]

列视角

若将系数矩阵按列切分，则每一列表示一个向量，方程组的解 (x, y) 表示每个列向量以 x, y 为权重的线性组合刚好形成 b 向量。这个就是方程组 \(Ax=b\) 有解的条件：b 在 A 的列空间，此时，x 为 列向量的组合系数。 \[ x \begin{bmatrix} \color{Green} 2 \\ \color{Green} -1 \end{bmatrix} + y \begin{bmatrix} \color{Turquoise} -1 \\ \color{Turquoise} 2 \end{bmatrix} = \begin{bmatrix} \color{Red} 0\\ \color{Red} 3 \end{bmatrix} \]

三元方程组

行视角

对于三元方程组行视角来说，每一行的方程组成一个三维空间中的一个平面。解是三个平面的交点，通常来说为一个点。 \[ \begin{bmatrix} \color{magenta} m_{11} & \color{magenta} m_{12} & \color{magenta} m_{13} \\ \color{cyan} m_{21} & \color{cyan} m_{22} & \color{cyan} m_{23} \\ \color{green} m_{31} & \color{green} m_{32} & \color{green} m_{33} \\ \end{bmatrix} \begin{bmatrix} w_{1} \\ w_{2} \\ w_{3} \end{bmatrix}= \begin{bmatrix} \color{magenta} y_{1} \\ \color{cyan} y_{2} \\ \color{green} y_{3} \end{bmatrix} \]

图片来自 Introduction to Linear Algebra for Applied Machine Learning with Python，解为一直线而非一个点。

列视角

三元列视角下，A 的列向量为三维平面的一个向量，x（下图为 w） 表示每个列向量取多少倍数可以组成 b 向量。 \[ \color{cyan} w_1 \color{black} \begin{bmatrix} \color{cyan} m_{11} \\ \color{cyan} m_{21} \\ \color{cyan} m_{31} \end{bmatrix}+ \color{magenta} w_2 \color{black} \begin{bmatrix} \color{magenta} m_{12} \\ \color{magenta} m_{22} \\ \color{magenta} m_{32} \end{bmatrix}+ \color{green} w_3 \color{black} \begin{bmatrix} \color{green} m_{13} \\ \color{green} m_{23} \\ \color{green} m_{33} \end{bmatrix}= \begin{bmatrix} \color{red} y_{1} \\ \color{red} y_{2} \\ \color{red} y_{3} \end{bmatrix} \]

解方程的步骤

总结了二元三元方程组的行和列视角后，我们回顾求解方程的具体步骤。用两个过程，消元和回代便可以解得方程。

• 消元的目的是将系数矩阵表示成变量依次依赖的上矩阵形式 (Upper Triangular Matrix)。 \[ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \\ \end{bmatrix} \underrightarrow{Eliminate} \begin{bmatrix} \color{red} a_{11} & a_{12} & \cdots & a_{1n} \\ \color{gray} 0 & \color{red} a'_{22} & \cdots & a'_{2n} \\ \color{gray} 0 & \vdots & \ddots & \vdots \\ \color{gray} 0 & \color{gray} 0 & \color{gray} 0 & \color{red} a'_{nn} \\ \end{bmatrix} \]

• 回代则在上矩阵的基础上依次解得每个分量的值。

三元方程示例

举个三元方程组为例 \[ \begin{cases} x&+2y&+z&=2 \\ 3x&+8y&+z&=12 \\ &+4y&+z&=2 \end{cases} \]

写成矩阵形式为

\[ \begin{bmatrix}1&2&1\\3&8&1\\0&4&1\end{bmatrix}\begin{bmatrix}x\\y\\z\end{bmatrix}=\begin{bmatrix}2\\12\\2\end{bmatrix} \]

消元过程

在消元过程中，有两类操作，一是将上一行乘以某系数后被下一行减去，依次消除这一行的元。第二类操作是交换当前行和后面某行，行交换用于当某行对应的元已经为0的情况下。

以上述三元方程为例，第二行消元的具体过程为第二行减去3倍的第一行。接着，再进行第三行消元。

\[ \left[\begin{array}{c|c}A&b\end{array}\right] = \left[\begin{array}{ccc|c} \color{red} 1&2&1&2 \\ 3&8&1&12 \\ 0&4&1&2 \end{array}\right] \xrightarrow{row_2-3row_1} \left[\begin{array}{ccc|c} \color{red} 1&2&1&2 \\ \color{blue} \boxed{0} & \color{red} 2& \color{blue} -2& \color{blue}6 \\ 0&4&1&2 \end{array}\right] \xrightarrow{row_3-2row_2} \left[\begin{array}{ccc|c} \color{red}1&2&1&2 \\ 0& \color{red} 2&-2&6 \\ \color{blue} 0 & \color{blue} \boxed{0} & \color{red} 5& \color{blue}-10 \end{array}\right] \]

最终，上矩阵为 \[ U=\begin{bmatrix} \color{red} 1 & 2& 1 \\ \color{gray} 0 & \color{red} 2 & -2 \\ \color{gray} 0 & \color{gray} 0& \color{red} 5 \end{bmatrix} \]

回代过程

回代过程比较直白，由上矩阵 U 对应方程组 \[ \begin{cases} \color{red} x&+2y&+z&=2 \\ & \color{red} +2y &-2z&=6 \\ && \color{red} +5z&=-10 \end{cases} \]

自下向上，容易解得 \[ \begin{cases} \color{red}x \color{black}= 2 \\ \color{red}y \color{black}= 1 \\ \color{red}z \color{black}= -2 \end{cases} \]

消元的行视角意义

注意到消元时的两类操作都不改变系数矩阵的行空间，只是改变了行空间的线性组合方式。

由于每一行代表一个拘束子空间，因此每一次消元改变了该行的拘束子空间。

举个例子，对于二元方程组和行空间 \[ \begin{bmatrix} \color{blue} 2 x & \color{blue} -1 y \\ \color{green} -1 x & \color{green} 2 y \end{bmatrix} = \begin{bmatrix} \color{blue} 0 \\ \color{green} 3 \end{bmatrix} \]

对应了两条直线

第二行的 x 消除后其几何意义为：蓝色直线不变，绿色直线从包含 x 的成分变成不含 x 成分，并且维持交点 （1, 2）不变。 \[ \begin{bmatrix} \color{blue} 2 x & \color{blue} -1 y \\ & \color{green} {3 \over 2} y \end{bmatrix} = \begin{bmatrix} \color{blue} 0 \\ \color{green} 3 \end{bmatrix} \]

最后还有一个问题可以思考一下：消元对于列视角的几何意义是什么呢？

在 MIT Gilbert Strang 教授经典的线性代数课程第三节中详细解释了矩阵相乘的5种方法及其理解，本文用彩色的 Latex 公式呈现出来，便于大家回顾复习。

01 - 矩阵乘法的五种理解

02 - 解方程组的意义和过程

假定 \(A \times B =C\) 中， A 是 \(m\) 行 \(n\) 列的矩阵， \(B\) 是 \(n\) 行 \(p\) 列的矩阵，\(C\) 为 \(m\) 行 \(p\) 列矩阵：

\[ A_{m \cdot n} \times B_{n \cdot p} = C_{m \cdot p} \]

B矩阵的列组合：列列切分

这是最经典的理解方式，沿袭了第一部分 \(A \vec x= \vec b\) 的方式。

回顾可以将 \(A \vec x\) 视为 \(A\) 的列向量关于每个 \(x\) 分量的线性组合。

那么 \(A B\) 相乘可以理解为将矩阵 \(B\) 按列切分成列向量，即 \[ B = [\vec B_{col_1} \vec B_{col_2} \cdots \vec B_{col_p} ] \]

如此，结果矩阵的第 \(j\) 列就是 \(A \vec B_{col_j}\)：\(A\) 的列向量关于 $ B_{col_j} $ 的线性组合。

由于我们将 \(A\) 和 \(B\) 都按列来切，这种方式可以助记成 列列 切分。

A矩阵行组合：行行切分

同样的，对应与 \(A\) 右乘向量等同于\(A\)列的组合，\(A\) 左乘行向量等同于 \(A\) 行的组合： \[ [x_1, x_2, ..., x_n]\begin{bmatrix} A_{row_1} \\ A_{row_2} \\ \vdots \\A_{row_n} \end{bmatrix} = x_1 \cdot A_{row_1} + x_2 \cdot A_{row_2} +...+ x_n \cdot A_{row_n} \] 其结果是一个行向量。

那么 \(A B\) 相乘可以理解为将矩阵 \(A\) 按行切分成行向量，即 \[ A = \begin{bmatrix} A_{row_1} \\ A_{row_2}\\ \vdots \\ A_{row_m} \end{bmatrix} \]

如此，结果矩阵的第 \(i\) 行就是 \(\vec A_{row_i} B\)：\(B\) 的行向量关于 $ A_{row_i} $ 的线性组合。

这种方式可以助记成 行行 切分。

A行 x B列 点乘：行列切分

如 \(A\) 矩阵按行切，\(B\) 矩阵按列切，即住记成 行列 切分时如下所示。

\[ \begin{eqnarray*} && A_{m \cdot n} \times B_{n \cdot p} \\ &=& \begin{bmatrix} A_{row_1} \\ A_{row_2} \\ \vdots \\ A_{row_i} \\ \vdots \\ \end{bmatrix}_{m \cdot n} \begin{bmatrix} B_{col_1} & B_{col_2} & \cdots & B_{col_p} \end{bmatrix}_{n \cdot p} \\ &=& \begin{bmatrix} &\vdots&\\& A_{row_i}&\\&\vdots& \end{bmatrix}_{m \cdot n} \begin{bmatrix} &&\\\cdots& B_{col_j} & \cdots\\&& \end{bmatrix}_{n \cdot p} \\ &=& \begin{bmatrix} &\vdots&\\\cdots&c_{ij}&\cdots\\& \vdots \end{bmatrix}_{m \cdot p} \end{eqnarray*} \]

行乘以列即列向量点乘，结果是一个标量。因此 \(c_{ij}\) 为结果矩阵 \(C\) 的第 \(i\) 行 \(j\) 列的值。

\[ c_{ij}=A_{row_i} \cdot B_{col_j} = \begin{bmatrix} a_{i1} & a_{i2} & \cdots & a_{in} \end{bmatrix} \begin{bmatrix} b_{1j} \\ b_{2j} \\ \vdots \\ b_{nj} \end{bmatrix} =\sum_{k=i}^na_{ik}b_{kj} \]

A列 x B行 矩阵和：列行切分

最后，也可也按列行来切分。注意列乘以行时的结果是一个矩阵。

分块相乘

第五种方式是分块相乘，可以认为是点乘理解下的扩展。

快速幂运算是一种利用位运算和DP思想求的\(x^n\)的数值算法，它将时间复杂度\(O(n)\)降到\(O(log(n))\)。快速幂运算结合矩阵乘法，可以巧解不少DP问题。本篇会由浅入深，从最基本的快速幂运算算法，到应用矩阵快速幂运算解DP问题，结合三道Leetcode题目来具体讲解。

Leetcode 50. Pow(x, n) (Medium)

Leetcode 50. Pow(x, n) 是实数的快速幂运算问题，题目如下。

Implement pow(x, n), which calculates x raised to the power n (i.e. \(x^n\)).

Example 1:

1
2

Input: x = 2.00000, n = 10
Output: 1024.00000

Example 2:

1
2

Input: x = 2.10000, n = 3
Output: 9.26100

Example 3:

1
2
3

Input: x = 2.00000, n = -2
Output: 0.25000
Explanation: 2-2 = 1/22 = 1/4 = 0.25

快速幂运算解法分析

假设n是32位的int类型，将n写成二进制形式，那么n可以写成最多32个某位为 1（第k位为1则值为\(2^k\)）的和。那么\(x^n\)最多可以由32个 \(x^{2^k}\)的乘积组合，例如：

快速幂运算的特点就是通过32次循环，每次循环根据上轮\(x^{2^k}\)的值进行平方后得出这一轮的值：\(x^{2^k} \times x^{2^k} = x^{2^{k+1}}\)，即循环计算出如下数列

在循环时，如果n的二进制形式在本轮对应的位的值是1，则将这次结果累乘计入最终结果。

下面是python 3 的代码，由于循环为32次，所以容易看出算法复杂度为 \(O(log(n))\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# AC
# Runtime: 32 ms, faster than 54.28% of Python3 online submissions for Pow(x, n).
# Memory Usage: 14.2 MB, less than 5.04% of Python3 online submissions for Pow(x, n).

class Solution:
    def myPow(self, x: float, n: int) -> float:
        ret = 1.0
        i = abs(n)
        while i != 0:
            if i & 1:
                ret *= x
            x *= x
            i = i >> 1
        return 1.0 / ret if n < 0 else ret

对应的 Java 的代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

// AC
// Runtime: 1 ms, faster than 42.98% of Java online submissions for Pow(x, n).
// Memory Usage: 38.7 MB, less than 48.31% of Java online submissions for Pow(x, n).

class Solution {
    public double myPow(double x, int n) {
        double ret = 1.0;
        long i = Math.abs((long) n);
        while (i != 0) {
            if ((i & 1) > 0) {
                ret *= x;
            }
            x *= x;
            i = i >> 1;
        }

        return n < 0 ? 1.0 / ret : ret;
    }
}

矩阵快速幂运算

快速幂运算也可以应用到计算矩阵的幂，即上面的x从实数变为方形矩阵。实现上，矩阵的幂需要矩阵乘法：$ A_{r c} B_{c p}$ ，Python中可以用numpy的 np.matmul(A, B)来完成，而Java版本中我们手动实现简单的矩阵相乘算法，从三重循环看出其算法复杂度为\(O(r \times c \times p)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

public int[][] matrixProd(int[][] A, int[][] B) {
    int R = A.length;
    int C = B[0].length;
    int P = A[0].length;
    int[][] ret = new int[R][C];
    for (int r = 0; r < R; r++) {
        for (int c = 0; c < C; c++) {
            for (int p = 0; p < P; p++) {
                ret[r][c] += A[r][p] * B[p][c];
            }
        }
    }
    return ret;
}

Leetcode 509. Fibonacci Number (Easy)

有了快速矩阵幂运算，我们来看看如何具体解题。Fibonacci问题作为最基本的DP问题，在上一篇Leetcode 679 24 Game 的 Python 函数式实现中我们用python独有的yield来巧解，这次再拿它来做演示。

The Fibonacci numbers, commonly denoted F(n) form a sequence, called the Fibonacci sequence, such that each number is the sum of the two preceding ones, starting from 0 and 1. That is,

1
2

F(0) = 0,   F(1) = 1
F(N) = F(N - 1) + F(N - 2), for N > 1.

Given N, calculate F(N).

Example 1:

1
2
3

Input: 2
Output: 1
Explanation: F(2) = F(1) + F(0) = 1 + 0 = 1.

Example 2:

1
2
3

Input: 3
Output: 2
Explanation: F(3) = F(2) + F(1) = 1 + 1 = 2.

Example 3:

1
2
3

Input: 4
Output: 3
Explanation: F(4) = F(3) + F(2) = 2 + 1 = 3.

转换为矩阵幂运算

AC代码

有了上面的矩阵幂公式，代码稍作改动即可。Java 版本代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

/**
 * AC
 * Runtime: 0 ms, faster than 100.00% of Java online submissions for Fibonacci Number.
 * Memory Usage: 37.9 MB, less than 18.62% of Java online submissions for Fibonacci Number.
 *
 * Method: Matrix Fast Power Exponentiation
 * Time Complexity: O(log(N))
 **/
class Solution {
    public int fib(int N) {
        if (N <= 1) {
            return N;
        }
        int[][] M = {{1, 1}, {1, 0}};
        // powers = M^(N-1)
        N--;
        int[][] powerDouble = M;
        int[][] powers = {{1, 0}, {0, 1}};
        while (N > 0) {
            if (N % 2 == 1) {
                powers = matrixProd(powers, powerDouble);
            }
            powerDouble = matrixProd(powerDouble, powerDouble);
            N = N / 2;
        }

        return powers[0][0];
    }

    public int[][] matrixProd(int[][] A, int[][] B) {
        int R = A.length;
        int C = B[0].length;
        int P = A[0].length;
        int[][] ret = new int[R][C];
        for (int r = 0; r < R; r++) {
            for (int c = 0; c < C; c++) {
                for (int p = 0; p < P; p++) {
                    ret[r][c] += A[r][p] * B[p][c];
                }
            }
        }
        return ret;
    }

}

Python 3的numpy.matmul() 版本代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# AC
# Runtime: 256 ms, faster than 26.21% of Python3 online submissions for Fibonacci Number.
# Memory Usage: 29.4 MB, less than 5.25% of Python3 online submissions for Fibonacci Number.

class Solution:

    def fib(self, N: int) -> int:
        if N <= 1:
            return N

        import numpy as np
        F = np.array([[1, 1], [1, 0]])

        N -= 1
        powerDouble = F
        powers = np.array([[1, 0], [0, 1]])
        while N > 0:
            if N % 2 == 1:
                powers = np.matmul(powers, powerDouble)
            powerDouble = np.matmul(powerDouble, powerDouble)
            N = N // 2

        return powers[0][0]

或者也可以直接调用numpy.matrix_power() 代替手动的快速矩阵幂运算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# AC
# Runtime: 116 ms, faster than 26.25% of Python3 online submissions for Fibonacci Number.
# Memory Usage: 29.2 MB, less than 5.25% of Python3 online submissions for Fibonacci Number.

class Solution:

    def fib(self, N: int) -> int:
        if N <= 1:
            return N

        from numpy.linalg import matrix_power
        import numpy as np
        F = np.array([[1, 1], [1, 0]])
        F = matrix_power(F, N - 1)

        return F[0][0]

Leetcode 1411. Number of Ways to Paint N × 3 Grid (Hard)

下面来看一道稍难一点的DP问题，1411. Number of Ways to Paint N × 3 Grid。

You have a grid of size n x 3 and you want to paint each cell of the grid with exactly one of the three colours: Red, Yellow or Green while making sure that no two adjacent cells have the same colour (i.e no two cells that share vertical or horizontal sides have the same colour).

You are given n the number of rows of the grid.

Return the number of ways you can paint this grid. As the answer may grow large, the answer must be computed modulo 10^9 + 7.

Example 1:

1
2
3

Input: n = 1
Output: 12
Explanation: There are 12 possible way to paint the grid as shown:

Example 2:

1
2

Input: n = 2
Output: 54

Example 3:

1
2

Input: n = 3
Output: 246

Example 4:

1
2

Input: n = 7
Output: 106494

Example 5:

1
2

Input: n = 5000
Output: 30228214

标准DP解法

分析题目容易发现第i行的状态只取决于第i-1行的状态，第i行会有两种不同状态：三种颜色都有或者只有两种颜色。这个问题容易识别出是经典的双状态DP问题，那么我们定义dp2[i]为第i行只有两种颜色的数量，dp3[i]为第i行有三种颜色的数量。

先考虑dp3[i]和i-1行的关系。假设第i行包含3种颜色，即dp3[i]，假设具体颜色为红，绿，黄，若i-1行包含两种颜色（即dp2[i-1]），此时dp2[i-1]只有以下2种可能：

因此，dp3[i]= dp2[i-1] * 2 + dp3[i-1] * 2。

同理，若第i行包含两种颜色，即dp2[i]，假设具体颜色为绿，黄，绿，若i-1行是两种颜色（dp2[i-1]），此时dp2[i-1]有如下3种可能：

因此，dp2[i] = dp2[i-1] * 3 + dp3[i-1] * 2。 初始值dp2[1] = 6，dp3[1] = 6，最终答案为dp2[i] + dp3[i]。

很容易写出普通DP版本的Python 3代码，时间复杂度为\(O(n)\)。

1
2
3
4
5
6
7
8
9
10
11
12
13

# AC
# Runtime: 36 ms, faster than 98.88% of Python3 online submissions for Number of Ways to Paint N × 3 Grid.
# Memory Usage: 13.9 MB, less than 58.66% of Python3 online submissions for Number of Ways to Paint N × 3 Grid.

class Solution:
    def numOfWays(self, n: int) -> int:
        MOD = 10 ** 9 + 7
        dp2, dp3 = 6, 6
        n -= 1
        while n > 0:
            dp2, dp3 = (dp2 * 3 + dp3 * 2) % MOD, (dp2 * 2 + dp3 * 2) % MOD
            n -= 1
        return (dp2 + dp3) % MOD

快速矩阵幂运算解法

代码几乎和Fibonacci一模一样，仅仅多了mod 计算。下面是Java版本。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

/**
AC
Runtime: 0 ms, faster than 100.00% of Java online submissions for Number of Ways to Paint N × 3 Grid.
Memory Usage: 35.7 MB, less than 97.21% of Java online submissions for Number of Ways to Paint N × 3 Grid.
**/

class Solution {
    public int numOfWays(int n) {
        long MOD = (long) (1e9 + 7);
        long[][] ret = {{6, 6}};
        long[][] m = {{3, 2}, {2, 2}};
        n -= 1;
        while(n > 0) {
            if ((n & 1) > 0) {
                ret = matrixProd(ret, m, MOD);
            }
            m = matrixProd(m, m, MOD);
            n >>= 1;
        }
        return (int) ((ret[0][0] + ret[0][1]) % MOD);

    }

    public long[][] matrixProd(long[][] A, long[][] B, long MOD) {
        int R = A.length;
        int C = B[0].length;
        int P = A[0].length;
        long[][] ret = new long[R][C];
        for (int r = 0; r < R; r++) {
            for (int c = 0; c < C; c++) {
                for (int p = 0; p < P; p++) {
                    ret[r][c] += A[r][p] * B[p][c];
                    ret[r][c] = ret[r][c] % MOD;
                }
            }
        }
        return ret;
    }

}

Python 3实现为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# AC
# Runtime: 88 ms, faster than 39.07% of Python3 online submissions for Number of Ways to Paint N × 3 Grid.
# Memory Usage: 30.2 MB, less than 11.59% of Python3 online submissions for Number of Ways to Paint N × 3 Grid.

class Solution:
    def numOfWays(self, n: int) -> int:
        import numpy as np

        MOD = int(1e9 + 7)
        ret = np.array([[6, 6]])
        m = np.array([[3, 2], [2, 2]])

        n -= 1
        while n > 0:
            if n % 2 == 1:
                ret = np.matmul(ret, m) % MOD
            m = np.matmul(m, m) % MOD
            n = n // 2
        return int((ret[0][0] + ret[0][1]) % MOD)