最近重新学习了一下线性代数,读完了《程序猿的数学3 线性代数》,为深刻理解内容,打算实现一个矩阵库,并写Python模块。计划涉及如下部分:

  • 向量内积
  • 矩阵乘法
  • LU分解
  • 线性方程组求解
  • 行列式计算
  • 逆矩阵
  • 实对称矩阵Jocobi方法(求解特征值)
  • QR分解

目前项目地址:https://github.com/netcan/LinAlg

关于Python模块的实现,可以参考我的这一篇文章:从头开始实现一个线性代数库:Python模块篇

实现

内部统一采用double数据类型。

LU分解

LU分解将矩阵A分解为两个矩阵的乘积$A=LU$,其中$L$为下三角矩阵,而$U$为上三角矩阵,$LU$的形式如下,实现的时候采用一个矩阵来存储$LU$的结果。
$$
L = \left(
\begin{matrix}
1 & 0 & 0\\
x & 1 & 0\\
x & x & 1\end{matrix}
\right)
U = \left(
\begin{matrix}
x & x & x & x\\
0 & x & x & x\\
0 & 0 & x & x\end{matrix}
\right)
$$

LU分解可用来方便的求解线性方程组、行列式的值以及逆矩阵。

分解方法:

将L, U改写成分块矩阵的形式:
$$
L = (l_1, \cdots, l_s), U = \left(\begin{matrix}u_1^T\\\vdots\\u_s^T\end{matrix}\right)
$$

$$
A = l_1u_1^t + \cdots + l_su_s^T
$$

对$A\equiv A(1)=(a_{ij}(1))$,令
$$
l_1 = \dfrac{1}{a_{11}(1)}\left(\begin{matrix}a_{11}(1)\\\vdots\\a_{m1}(1)\end{matrix}\right) \quad
u_1^t = (a_{11}(1), \cdots, a_{1n}(1))
$$


$$
A(1) -l_1u_1^T =
\left(
\begin{array}{c|ccc}
0 & 0 & \cdots & 0 \\ \hline
0 & & & \\
\vdots & & A(2) & \\
0 & & & \\
\end{array}
\right)
$$

$$
A(1) = l_1u_1^T +
\left(
\begin{array}{c|ccc}
0 & 0 & \cdots & 0 \\ \hline
0 & & & \\
\vdots & & A(2) & \\
0 & & & \\
\end{array}
\right)
$$

对$A(s), s = 2, 3,\cdots$进行和$A(1)$一样的操作,最终得到零矩阵分解完毕可以得到$L$, $U$。

核心代码如下:

Matrix Matrix::LUdecomp() const {
    Matrix LU(*this); // 最终结果存放到一个矩阵中
    size_t s = std::min(getRowSize(), getColSize());
    for(size_t k = 0; k < s; ++k) {
        double x = 1.0 / LU.data[k][k];
        // l1 = 1/a_{11} * ...
        for(size_t i = k + 1; i < getRowSize(); ++i)
            LU.data[i][k] *= x;
        // A -= lu^t,
        for(size_t i = k + 1; i < getRowSize(); ++i)
            for(size_t j = k + 1; j < getColSize(); ++j)
                LU.data[i][j] -= LU.data[i][k] * LU.data[k][j];
    }
    return LU;
}

若分解过程中出现除数为零,这需要添加一个置换矩阵$P$,分解得到$A = PLU$的形式,详见Matrix::PLUdecomp()

行列式计算

对$A$进行$A=PLU$分解后,有
$$
det(A) = det(PLU) = det(P)det(L)det(U) = det(P)det(U)
$$

值即为$U$对角线上元素的乘积。 核心代码如下:

double Matrix::det() const {
    assert(getRowSize() == getColSize());
    auto PLU = PLUdecomp();
    const auto &P = PLU.first;
    const auto &LU = PLU.second;

    bool sign = false; // 负号
    for(size_t i = 0; i < PLU.first.size(); ++i)
        for(size_t j = i + 1; j < PLU.first.size(); ++j)
            if(PLU.first[i] > PLU.first[j]) sign = !sign;
    double ret = sign?-1.0:1.0;
    // det(A) = det(PLU) = det(P)det(L)det(U) = det(P)det(U)
    for(size_t i = 0; i < LU.data.size(); ++i)
        ret *= LU.data[P[i]][i];
    return ret;
}

求解线性方程组$Ax=y$

对$A$进行$A=PLU$分解后,有
$$
Ax = LUx = L(Ux) = Lz = y\\
Lz = \left(
\begin{matrix}
1 & & & \\
a_{21} & 1 & & \\
a_{31} & a_{32} & 1 & \\
a_{41} & a_{42} & a_{43} & 1\\
\end{matrix}
\right)
\left(\begin{matrix}z_1\\z_2\\z_3\\z_4\end{matrix}\right)
=
\left(\begin{matrix}y_1\\y_2\\y_3\\y_4\end{matrix}\right) = y\\
Ux = \left(
\begin{matrix}
a_{11} & a_{12} & a_{13} & a_{14} \\
& a_{22} & a_{23} & a_{24}\\
& & a_{33} & a_{34} \\
& & &a_{44}\\
\end{matrix}
\right)
\left(\begin{matrix}x_1\\x_2\\x_3\\x_4\end{matrix}\right)
=
\left(\begin{matrix}z_1\\z_2\\z_3\\z_4\end{matrix}\right) = z\\
$$

先用$y$解出$z$,在用$z$解出$x$。在解$Lz=y$时,从头开始计算;而在解$Ux = z$时,从尾开始计算。 核心代码如下:

Vector Matrix::LUsolve_L(const Matrix &LU, const Vector &y) const {
    // 求出Lz = y的z
    Vector z(y.getSize());
    for(size_t i = 0; i < z.getSize(); ++i) {
        z.data[i] = y.data[i];
        for (size_t j = 0; j < i; ++j)
            z.data[i] -= LU.data[i][j] * z.data[j];
    }
    return z;
}

Vector Matrix::LUsolve_U(const Matrix &LU, const Vector &z) const {
    // 求出Ux = z的x
    Vector x(z.getSize());
    for(int i = z.getSize() - 1; i >=0; --i) {
        x.data[i] = z.data[i];
        for(int j = i + 1; j < z.getSize(); ++j)
            x.data[i] -= LU.data[i][j] * x.data[j];
        x.data[i] /= LU.data[i][i];
    }
    return x;
}

Vector Matrix::LUsolve(const Vector &y) const {
    // Ax = y
    assert(y.getSize() == getRowSize());
    assert(! isZero(det()) && "det is equal 0");
    Matrix LU = std::move(LUdecomp());
    Vector x(y.getSize());
    x = std::move(LUsolve_L(LU, y));
    x = std::move(LUsolve_U(LU, x));
    return x;
}

求解逆矩阵

对于$AX=AA^{-1}=I$,有
$$
A(x_1, \cdots, x_n) = (e_1, \cdots, e_n)展开得\\
Ax_1 = e_1, \cdots, Ax_n = e_n
$$

利用前面求解线性方程组的方法,得到$X=A^{-1}$。核心代码如下:

Matrix Matrix::inv() const {
    // A(x_1, ..., x_n) = (e_1, ..., e_n)
    assert(! isZero(det()) && "matrix is not invertible");
    size_t n = getRowSize();
    Matrix LU = std::move(LUdecomp());
    Matrix ret(n, n);
    Vector e(n), x(n);
    e[0] = 1.0;
    for(size_t j = 0; j < n; ++j) {
        if(j >= 1) std::swap(e[j], e[j-1]);
        // 解出Ax = e
        x = std::move(LUsolve_L(LU, e));
        x = std::move(LUsolve_U(LU, x));
        for(size_t i = 0; i < n; ++i)
            ret.data[i][j] = x[i];
    }
    return ret;
}

利用Jacobi算法求解实对称矩阵特征值

对于实对称矩阵$A$,将其对角化,有$P^{-1}AP=D$,这里的$D$为对角矩阵,对角线上的值即为特征值,$P$为特征向量组成的矩阵。

而解实对角矩阵的特征值算法有Jacobi,主要通过不断平面旋转变化,使得矩阵$A$朝着$D$的方向变换。

定义如下平面旋转矩阵(正交矩阵)$R(\theta, p, q)$:
$$
R(\theta, p, q) =
\begin{array}{cc}
& \begin{array}{ccccccccccc}
& & & 第p列 & & & & 第q列 & & & &
\end{array}\\
\begin{array}{c}
\\
\\
\\
第p行\\
\\
\\
\\
第q行\\
\\
\\
\\
\end{array}
& \left(
\begin{array}{ccccccccccc}
1 & & & & & & & & & & \\
& \ddots & & & & & & & & & \\
& & 1 & & & & & & & & \\
& & &\cos\theta & & & -\sin\theta & & & & \\
& & & & 1 & & & & & & \\
& & & & &\ddots & & & & & \\
& & & & & & 1 & & & & \\
& & & \sin\theta & & & \cos\theta & & & & \\
& & & & & & & & 1 & & \\
& & & & & & & & & \ddots & \\
& & & & & & & & & & 1 \\
\end{array}
\right) \\
\end{array}
$$

易知$R^T R = I$,所以有$R^{-1} = R^T$。通过正交矩阵将对称矩阵$A$不断旋转变换进行对角化,$D = R(\theta, p, q)^{-1}AR(\theta, p, q) = R(\theta, p, q)^TAR(\theta, p, q)$。对$A$左乘$R^T$后,$p, q$两行将发生变化,而右乘$R$后,$p,q$两列发生变化。

$$
R(\theta, p, q)^T
\left(\begin{matrix}a_{1j}\\\vdots\\a_{pj}\\\vdots\\a_{qj}\\\vdots\\a_{nj}\end{matrix}\right)
= \left(\begin{matrix}a_{1j}\\\vdots\\a_{pj}\cos\theta + a_{qj}\sin\theta \\\vdots\\ -a_{pj}\sin\theta + a_{qj}\cos\theta \\\vdots\\a_{nj}\end{matrix}\right) =
\left(\begin{matrix}a_{1j}\\\vdots\\a_{pj}’\\\vdots\\a_{qj}’\\\vdots\\a_{nj}\end{matrix}\right) \\
\begin{array}{ccc}
(a_{i1}, \cdots, a_{ip}, \cdots, a_{iq}, \cdots, a_{in}) R(\theta, p, q) &=& (a_{i1},& \cdots,& a_{ip} \cos\theta + a_{iq} \sin\theta,& \cdots,& -a_{ip}\sin\theta + a_{iq}\cos\theta,& \cdots,& a_{in})\\
& & (a_{i1},& \cdots,& a_{ip}’,& \cdots,& a_{iq}’,& \cdots,& a_{in})
\end{array}
$$

将p, q两行两列的将四个交叉点的结果写出来,由于只有这四个交叉点同时受到行列影响,所以它们的公式略有不同(由于$A$是对称矩阵有$a_{pq} = a_{qp}$):
$$
\begin{array}{ccl}
a_{pp}’ &=& (a_{pp}\cos\theta + a_{qp}\sin\theta)\cos\theta + (a_{pq}cos\theta + a_{qq}\sin\theta)\sin\theta\\
& = & a_{pp}\cos^2\theta + a_{qq}\sin^2\theta + 2a_{pq}\sin\theta\cos\theta\\
a_{pq}’ &=& -(a_{pp}\cos\theta + a_{qp}\sin\theta)\sin\theta + (a_{pq}cos\theta + a_{qq}\sin\theta)\cos\theta\\
& = & a_{pq}(\cos^2\theta - \sin^2\theta) + (a_{qq} - a_{pp})\sin\theta\cos\theta\\
a_{qp}’ &=& (-a_{pp}\cos\theta + a_{qp}\sin\theta)\cos\theta + (-a_{pq}cos\theta + a_{qq}\sin\theta)\sin\theta\\
& = & a_{pq}(\cos^2\theta - \sin^2\theta) + (a_{qq} - a_{pp})\sin\theta\cos\theta\\
a_{qq}’ &=& -(-a_{pp}\cos\theta + a_{qp}\sin\theta)\sin\theta + (-a_{pq}cos\theta + a_{qq}\sin\theta)\cos\theta\\
& = & a_{pp}\sin^2\theta + a_{qq}\cos^2\theta - 2a_{pq}\sin\theta\cos\theta\\
\end{array}
$$

经过反复旋转变换后,矩阵$A$依旧保持着对称性$a_{pq} = a_{qp}$,依次选择$a_{pq}$,令$a_{pq}$旋转变化为$0$得到$\theta$,在这个$\theta$作用下更新其他$p, q$的行、列值。简而言之,每次将两行两列相交的四个交叉点中的斜对角交叉点变为0,最终趋向于对角矩阵。
$$
a_{pq}’ = a_{qp}’ = a_{pq}(\cos^2\theta - \sin^2\theta) + (a_{qq} - a_{pp})\sin\theta\cos\theta = 0\\
\Rightarrow \dfrac{a_{pq}}{a_{pp} - a_{qq}} = \dfrac{\sin\theta\cos\theta}{\cos^2\theta - \sin^2\theta} = \dfrac{1}{2}\tan2\theta
$$

通过$\tan2\theta$得到$\cos\theta$的值,有:
$$
1+\tan^{2}2\theta = \dfrac{\cos^{2}2\theta + \sin^{2}2\theta}{\cos^{2}2\theta} = \dfrac{1}{\cos^{2}{2\theta}}\\
\Rightarrow \cos2\theta = \dfrac{1}{\sqrt{1+\tan^{2}{2\theta}}} = \cos^2\theta - \sin^2\theta = 2\cos^2\theta - 1(0\leq \theta \leq \dfrac{\pi}{4})\\
\Rightarrow \cos\theta = \sqrt{\dfrac{1}{2}(1+\dfrac{1}{\sqrt{1+\tan^{2}2\theta}})}\\
\Rightarrow \sin\theta = \sqrt{1-\cos^2\theta}
$$

设$g(A) = \sum_{i}^na_{ii}^2$,当$g(A)$不在增大时,可视为$A$近似于对角矩阵$D$。核心代码如下:

void Matrix::R(double cosphi, size_t p, size_t q, bool T) {
    // 通过平面旋转进行相似变换
    double sinphi = sqrt(1 - cosphi * cosphi);
    size_t n = getRowSize();
    for(size_t k = 0; k < n; ++k) {
        double  ap =  data[T?p:k][T?k:p] * cosphi + data[T?q:k][T?k:q] * sinphi,
                aq = -data[T?p:k][T?k:p] * sinphi + data[T?q:k][T?k:q] * cosphi;
        data[T?p:k][T?k:p] = ap; data[T?q:k][T?k:q] = aq;
    }
    return;
}
vector<double> Matrix::Jacobi() const {
    assert(getRowSize() == getColSize());
    size_t n = getRowSize();
    // 判断是否对称
    for(size_t i = 0; i < n; ++i)
        for(size_t j = 0; j < n; ++j)
            assert(Eq(data[i][j], data[j][i]));
    vector<double> eigenvals(n);
    Matrix D(*this); // 对角化
    // g(A) = \sum a[i][i]^2
    auto g = [](const Matrix &A)->double {
        double ret = 0.0;
        for(size_t i = 0; i < A.getRowSize(); ++i)
            ret += A.data[i][i] * A.data[i][i];
        return ret;
    };
    // f(A) -> 0表明近似于对角化矩阵
    auto f = [](const Matrix &A)->double {
        double ret = 0.0;
        for(size_t i = 0; i < A.getRowSize(); ++i)
            for(size_t j = 0; j < A.getRowSize(); ++j)
                if(i!=j) ret += A.data[i][j] * A.data[i][j];
        return ret;
    };

    double cur_g = g(D), last_g = 0.0;

    while (! Eq(cur_g, last_g)) { // 当g(D)不变时,对角化完成
        for (size_t p = 0; p < n; ++p) { // 逐个处理a_{pq}
            for (size_t q = p + 1; q < n; ++q) {
                if (isZero(D.data[p][q])) continue;
                double  tan2phi = 2 * D.data[p][q] / (D.data[p][p] - D.data[q][q]),
                        cosphi = sqrt(0.5 * (1.0 + 1.0/sqrt(1+tan2phi * tan2phi)));
                // D = R^T A R,旋转变换
                D.R(cosphi, p, q, true);
                D.R(cosphi, p, q, false);
            }
        }
        last_g = cur_g;
        cur_g = g(D);
    }
    for(size_t i = 0; i < n; ++i) eigenvals[i] = D.data[i][i];
    return eigenvals;
}