动机、参考资料、涉及内容

动机

  • GPTQ 算法中涉及到一些关于矩阵的逆, 矩阵分解的内容

参考资料

待定

涉及内容

  • 线性代数定理(重结论, 适当证明)
  • numpy/torch 相关函数的使用

不涉及内容

待定

高斯消元法 (Gaussian elimination)

高斯消元法 (Gaussian elimination) 中的重要步骤是:

将一个矩阵的第 $i$ 行乘上一个系数 $k$ 加到第 $j$ 行上

\[\begin{bmatrix} a_{1,1} &a_{1,2} &\cdots &a_{1,n-1} &a_{1,n}\\ a_{2,1} &a_{2,2} &\cdots &a_{2,n-1} &a_{2,n}\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ a_{j,1} &a_{j,2} &\cdots &a_{j,n-1} &a_{j,n}\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ a_{m-1,1} &a_{m-1,2} &\cdots &a_{m-1,n-1} &a_{m-1,n}\\ a_{m,1} &a_{m,2} &\cdots &a_{m,n-1} &a_{m,n}\\ \end{bmatrix} \to \begin{bmatrix} a_{1,1} &a_{1,2} &\cdots &a_{1,n-1} &a_{1,n}\\ a_{2,1} &a_{2,2} &\cdots &a_{2,n-1} &a_{2,n}\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ k\times a_{i, 1}+a_{j,1} &k\times a_{i, 2}+a_{j,2} &\cdots &k\times a_{i, n-1}+a_{j,n-1} &k\times a_{i, n}+a_{j,n}\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ a_{m-1,1} &a_{m-1,2} &\cdots &a_{m-1,n-1} &a_{m-1,n}\\ a_{m,1} &a_{m,2} &\cdots &a_{m,n-1} &a_{m,n}\\ \end{bmatrix}\]

这个操作用矩阵乘法来表达可以写作

\[(\mathbf{I}_{(n,n)}+k\mathbf{E}_{j,i})A= \begin{bmatrix} 1 & 0 & \cdots &0 &0 \\ 0 & 1 & \cdots &0 &0 \\ \vdots &\vdots &\ddots &\vdots &\vdots \\ 0 & \dots & (j,i)=k, \dots &(j, j)=1 &0\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ 0 & 0 & \cdots &1 &0 \\ 0 & 0 & \cdots &0 &1 \\ \end{bmatrix} \begin{bmatrix} a_{1,1} &a_{1,2} &\cdots &a_{1,n-1} &a_{1,n}\\ a_{2,1} &a_{2,2} &\cdots &a_{2,n-1} &a_{2,n}\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ a_{j,1} &a_{j,2} &\cdots &a_{j,n-1} &a_{j,n}\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ a_{n-1,1} &a_{n-1,2} &\cdots &a_{n-1,n-1} &a_{n-1,n}\\ a_{m,1} &a_{m,2} &\cdots &a_{m,n-1} &a_{m,n}\\ \end{bmatrix}\]

其中 $\mathbf{E}_{j,i}$ 为 $(j, i)$ 元为 $1$, 其余元素均为 $0$ 的 $m\times m$ 方阵, 同样的道理,

将一个矩阵的第 $i$ 列乘上一个系数 $k$ 加到第 $j$ 列上, 用矩阵乘法表示是:

\[A(\mathbf{I}_{(m,m)}+k\mathbf{E}_{i, j})\]

其中 $\mathbf{E}_{i,j}$ 为 $(i, j)$ 元为 $1$, 其余元素均为 $0$ 的 $n\times n$ 方阵

小结论

(1)

\[\mathbf{E}_{j, i}A\mathbf{E}_{i, j}=a_{i,i}\mathbf{E}_{j, j}\]

(2) 假设 $A$ 是方阵, 且对角线元素 $a_{i, i}\neq 0$, 那么如下矩阵的第 $i$ 行与第 $i$ 列均为 $0$

\[A-\frac{1}{a_{i,i}}A_{:,i}A_{i,:}\]

注意:

\[A_{:,i}A_{i,:}=(A\mathbf{E}_{i, i})(\mathbf{E}_{i, i}A)=A\mathbf{E}_{i, i}A\]

高斯消元法最终是希望将(增广)矩阵转化为所谓的既约形行阶梯形矩阵, 除了上述变换外, 还需要使用到这两类操作:

将第 $i$ 行与第 $j$ 行交换

对应的矩阵变换为, 在矩阵 $A$ 的左边乘上一个矩阵:

\[\begin{align*} (\mathbf{I}-\mathbf{E}_{i,i}-\mathbf{E}_{j,j}+\mathbf{E}_{i,j}+\mathbf{E}_{j,i})A&= \begin{bmatrix} 1 &\cdots &\cdots &\cdots &\cdots &\cdots &0\\ \vdots &\ddots &\vdots &\ddots &\vdots &\ddots &\vdots\\ 0 & \cdots &(i,i)=0 &\cdots &(i, j)=1 &\cdots &0\\ \vdots &\ddots &\vdots &\ddots &\vdots &\ddots &\vdots\\ 0 & \cdots &(j, i)=1 &\cdots &(j, j)=0 &\cdots &0\\ \vdots &\ddots &\vdots &\ddots &\vdots &\ddots &\vdots\\ 0 &\cdots &\cdots &\cdots &\cdots &\cdots &1\\ \end{bmatrix}A \\ &= \begin{bmatrix} a_{1,1} &a_{1,2} &\cdots &a_{1,n-1} &a_{1,n}\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ {\color{red} a_{j,1} } &{\color{red} a_{j,2} }&\cdots & {\color{red} a_{j,n-1} }&{\color{red} a_{j,n} }\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ {\color{red} a_{i,1} } &{\color{red} a_{i,2} }&\cdots & {\color{red} a_{i,n-1} }&{\color{red} a_{i,n} }\\ \vdots &\vdots &\ddots &\vdots &\vdots \\ a_{m,1} &a_{m,2} &\cdots &a_{m,n-1} &a_{m,n}\\ \end{bmatrix} \end{align*}\]

将第 $i$ 行乘上一个数 $k$

\[(\mathbf{I}+(k-1)\mathbf{E}_{i,i})A\]

我们将这三类矩阵称为初等矩阵

\[\begin{matrix} P_{i,j}(k)=(\mathbf{I}+k\mathbf{E}_{i, j}) & \text{第}j\text{行乘以}k\text{加到第}i\text{行上} \\ P_{i,j}=(\mathbf{I}-\mathbf{E}_{i,i}-\mathbf{E}_{j,j}+\mathbf{E}_{i,j}+\mathbf{E}_{j,i})A& \text{交换第}i\text{行与第}j\text{行} \\ P_{i}(k) = (\mathbf{I}+(k-1)\mathbf{E}_{i,i}) &\text{第}i\text{行乘上}k \end{matrix}\]

行列式 (Determinant)

一个 $n$ 阶方阵的行列式定义为:

\[det(A)=\lvert A \rvert = \sum_{\sigma\in S_n}{sgn(\sigma)a_{1,\sigma(1)}\times\cdots\times a_{n, \sigma(n)}}\]

这里 $\sigma$ 是一个 ${1,\ldots,n}$ 的置换函数, $S_n$ 是所有的置换, 而 $sig(\sigma)$ 表示置换的符号, 如果 $\sigma$ 能通过偶数次交换变回恒等映射, 则符号为 $+1$, 否则为 $-1$. 从定义上看, 计算行列式的时间复杂度为 $n\times n!$

行列式的计算【待补充】

(1) 上三角或下三角矩阵的逆为对角线元素之积

(2) $det(AB)=det(A)det(B)$, 特别地, 因为高斯消元法对应的

矩阵的逆

定义矩阵 $A$ 的余子矩阵 (cofactor matrix) 为:

\[C=\begin{bmatrix} (-1)^{1+1}det(A_{-1,-1}), &\cdots, &(-1)^{1+n}det(A_{-1,-n}) \\ \vdots &\ddots &\vdots \\ (-1)^{n+1}det(A_{-n,-1}), &\cdots, &(-1)^{n+n}det(A_{-n,-n}) \\ \end{bmatrix}\]

其中 $A_{-i,-j}$ 表示除去 $A$ 的第 $i$ 行和第 $j$ 列构成的子矩阵.

定义矩阵 $A$ 的伴随矩阵 (adjoint matrix) 为余子矩阵的转置, 记作 $adj(A)$, 而矩阵的逆可以用伴随矩阵来表示:

\[A^{-1}=\frac{1}{det(A)}adj(A)=\frac{1}{det(A)}C^T\] 
import torch
from itertools import permutations

def sign(perm):
    # 使用逆序数计算符号
    s = 0
    n = len(perm)
    for i in range(n):
        for j in range(i+1, n):
            if perm[i] > perm[j]:
                s += 1
    return 1 if s % 2 == 0 else -1

def det(A):
    n = len(A)
    s = 0.
    for perm in permutations(range(n)):
        x = 1.
        for i, k in enumerate(perm):
            x *= A[i][k]
        s += x * sign(perm)
    return s

def submatrix(A, i, j):
    A1 = torch.concat([A[:i, :], A[i+1:, :]], axis=0)
    A1 = torch.concat([A1[:, :j], A1[:, j+1:]], axis=1)
    return A1

def cofactor(A):
    n = A.size(0)
    cof = torch.empty_like(A)
    for i in range(n):
        for j in range(n):
            sub = submatrix(A, i, j)
            s = -1 if (i+j) % 2 == 1 else 1
            cof[i][j] = s * det(sub)
    return cof

def adj(A):
    return cofactor(A).T

def inv(A):
    return 1 / det(A) * adj(A)

import torch
n = 4
A = torch.randn(n, n)

torch.allclose(torch.linalg.det(A), det(A))
torch.allclose(torch.linalg.inv(A), inv(A))

正定/半正定矩阵 (PSD/PD)

实正定矩阵对角线元素必为正数?

实半正定矩阵对角线元素非负? 分解后非零元的个数等于秩

实对称矩阵/正交/正定矩阵的分解形式

Cholesky 矩阵分解 (Cholesky decomposition)

PSD (Positive Semi-Definite matrix) 表示半正定矩阵, PD (Positive Definite matrix) 表示正定矩阵. 对于一个正定矩阵 $A$, 可以将 $A$ 进行如下矩阵分解:

\[A = LL^T,\quad L_{i,j}=0, \forall i < j\]

其中 $L$ 是一个下三角矩阵(对角线以上全为 $0$),

import torch
from torch.utils import benchmark

n = 256
k = 100

a = torch.randn((n, n))
M = a @ a.T + torch.eye(n)
y = torch.empty_like(M)
z = torch.empty_like(M)

t = benchmark.Timer(stmt="torch.inverse(M, out=y)", globals={"M": M, "y": y}).blocked_autorange(min_run_time=3).median
print(f"inverse ({n}, {n}) matrix, median time usage: {t*1000}ms")  # 0.71ms

t = benchmark.Timer(stmt="torch.cholesky_inverse(torch.linalg.cholesky(M), out=z)", globals={"M": M, "z": z}).blocked_autorange(min_run_time=3).median
print(f"cholesky_inverse ({n}, {n}) matrix, median time usage: {t*1000}ms")  # 0.60ms

Cholesky Algorithm

对于一个正定矩阵 $A$, 定义一系列 $n\times n$ 方阵: $A^{(1)},\ldots,A^{(n+1)}$, 其中 $A^{(1)}=A, A^{(n+1)}=\mathbf{I}$,

\[A^{(i)}=\begin{bmatrix} \mathbf{I}_{i-1} & 0 & \mathbf{0}^T \\ 0 & a_{i, i} & \mathbf{b}_i^T\\ \mathbf{0} &\mathbf{b}_i &B^{(i)} \end{bmatrix}\]

我们同样定义一系列下三角矩阵 $L^{(1)},\ldots L^{(n)}$

\[L^{(i)}=\begin{bmatrix} \mathbf{I}_{i-1} & 0 & \mathbf{0}^T \\ 0 & \sqrt{a_{i, i} } & \mathbf{0}^T\\ \mathbf{0} &\frac{1}{\sqrt{a_{i, i} } }\mathbf{b}_i & \mathbf{I}_{n-i} \end{bmatrix}\]

我们可以验证

\[A^{(i)}=L^{(i)}A^{(i+1)}(L^{(i)})^T\]

而

\[A^{(i+1)} = \begin{bmatrix} \mathbf{I}_{i-1} & 0 & \mathbf{0}^T \\ 0 & 1 & \mathbf{0}^T\\ \mathbf{0} &\mathbf{0} &B^{(i)}-\frac{1}{a_{i, i} }\mathbf{b}_i\mathbf{b}_i^T \end{bmatrix}= \begin{bmatrix} \mathbf{I}_{i} & 0 & \mathbf{0}^T \\ 0 & a_{i+1, i+1} & \mathbf{b}_{i+1}^T\\ \mathbf{0} &\mathbf{b}_{i+1} &B^{(i+1)} \end{bmatrix}\]

如此, 我们最终会得到 $A=L^{(1)}L^{(2)}\cdots L^{(n)}\mathbf{I}(L^{(n)})^T\cdots (L^{(2)})^T(L^{(1)})^T$

即定义 $L=L^{(1)}L^{(2)}\cdots L^{(n)}$, 我们有 $A=LL^T$, 也就是 Cholesky 分解形式.

并且我们可以验证上述分解过程与最终得到的下三角矩阵 $L$ 有如下直接联系:

\[L[i:, i] = L^{(i)}[i:, i],\quad i=1,...,n\]

备注: 以上式子中的 $a_{i,i}$ 不是原始矩阵 $A$ 的 $(i,i)$ 元($a_{1,1}$ 除外)