OI线性代数——矩阵

OI线性代数——矩阵

本章只讨论OI中的线性代数，过于高深的线性代数的知识请参考其他教材。

矩阵与向量

矩阵是一个数学概念，对应着计算机中的二维数组。例如：

$$A=\left[\begin{array}{ccc}1& 3& 5\\ 1& 2& 4\end{array}\right]$$

是一个$2\times 3$矩阵，由两行三列组成，通过$A[i,j]$来指定矩阵中第$i$行第$j$列的元素，例如$A[1,3]=3$。

有一类特殊的$n\times 1$矩阵叫做向量，例如：

$$V=\left[\begin{array}{c}1\\ 2\\ 3\end{array}\right]$$

是一个拥有$3$个元素的矩阵，此篇文章不再介绍计算几何中的向量。

方阵是一种$n\times n$的矩阵，例如：

$$B=\left[\begin{array}{cc}1& 3\\ 1& 2\end{array}\right]$$

是一个$2\times 2$的方阵，一般我们研究方阵的性质。

矩阵的运算

矩阵基本运算

• 矩阵加法：定义$A+B$，为两个相同大小的矩阵对应元素分别相加。
• 矩阵减法：定义$A-B$，为两个相同大小的矩阵对应元素分别相减。
• 矩阵数乘：定义$aA$，为矩阵$A$的每一个元素都乘以$a$。
• 矩阵转置：定义$A^T$，将每一个$A[i,j]$和$A[j,i]$互换位置。

矩阵乘法

定义$AB$，为$A$矩阵$a\times n$乘以矩阵$B$矩阵$n\times b$，结果矩阵为大小为$a\times b$。

完整定义：

$$AB[i,j]=\sum _{k=1}^{n}A[i,k]\times B[k,j]$$

矩阵乘法满足结合律，但不满足交换律。

定义单位矩阵为一个$n\times n$的方阵，其对角元素均为$1$，其他元素均为$0$。一个矩阵乘以单位矩阵不会改变这个矩阵。

矩阵的幂

定义矩阵的幂为$A^k=A\cdot A\dots A$，即$k$个$A$矩阵相乘。

矩阵相乘的朴素算法需要$O(k{n}^3)$时间，经过矩阵快速幂优化之后可以达到$O(n^3\log k)$的时间。

同理其他快速幂，我们也采用二进制快速幂的方式进行优化。

P3390 矩阵快速幂

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

#define FR freopen("in.txt", "r", stdin)
#define FW freopen("out.txt", "w", stdout)

int n;
ll k;

ll mod = 1000000007;

struct Matrix
{
    ll arr[105][105];

    Matrix()
    {
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= n; j++)
            {
                arr[i][j] = 0;
            }
        }
    }

    Matrix operator*(Matrix &o)
    {
        Matrix ans;
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= n; j++)
            {
                for (int k = 1; k <= n; k++)
                {
                    ans.arr[i][j] = (ans.arr[i][j] + (arr[i][k] * o.arr[k][j]) % mod) % mod;
                }
            }
        }
        return ans;
    }

    Matrix fpow(ll e)
    {
        Matrix temp = *this;
        Matrix ans;

        for (int i = 1; i <= n; i++)
            ans.arr[i][i] = 1;

        for (; e; e >>= 1, temp = temp * temp)
        {
            if (e & 1)
                ans = ans * temp;
        }

        return ans;
    }
};

int main()
{
    scanf("%d %lld", &n, &k);

    Matrix m;
    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            scanf("%lld", &m.arr[i][j]);
        }
    }

    m = m.fpow(k);

    for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= n; j++)
        {
            printf("%lld ", m.arr[i][j]);
        }
        printf("\n");
    }
    return 0;
}

线性递推

一个线性递推方程是一种数列，具有一些初始值$f(0),f(1),\dots ,f(k-1)$和一个递推式：

$$f(n)=c_{1}f(n-1)+c_{2}f(n-2)+\dots +c_{k}f(n-k)$$

这里$c_i$都是常数。普通的递推法的时间复杂度是$O(kn)$，而使用矩阵快速幂递推法的时间是$O(k^3\log n)$。

斐波那契数列递推

斐波那契数列有如下递推式：

$$\begin{gathered} f(0)=0 \\ f(1)=1 \\ f(n)=f(n-1)+f(n-2) \end{gathered}$$

我们定义数列向量为：

$$V_k=\left[\begin{array}{c}f(k)\\ f(k-1)\end{array}\right]$$

定义递推矩阵为：

$$V_{k+1}=T{V}_k$$

使用矩阵进行线性递推的过程需要构造出$T$矩阵。通过合理构造，我们有：

$$T=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$$

故有：

$$V_n=T^n\left[\begin{array}{c}1\\ 0\end{array}\right]$$

$T^n$可以在对数时间内计算出来。

更一般的递推

类比于上述定义，对于式子：

$$f(n)=c_{1}f(n-1)+c_{2}f(n-2)+\dots +c_{k}f(n-k)$$

我们可以构造出更一般的$T$矩阵为：

$$T=\left[\begin{array}{cccc}{c}_1& c_2& \cdots & c_k\\ 1& \cdots & 0& 0\\ ⋮& \ddots & ⋮& ⋮\\ 0& \cdots & 1& 0\end{array}\right]$$

P1939

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

#define FR freopen("in.txt", "r", stdin)
#define FW freopen("out.txt", "w", stdout)

int n;

ll mod = 1000000007;

struct Matrix
{
    ll arr[4][4];

    Matrix()
    {
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= n; j++)
            {
                arr[i][j] = 0;
            }
        }
    }

    Matrix operator*(const Matrix &o)
    {
        Matrix ans;
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= n; j++)
            {
                for (int k = 1; k <= n; k++)
                {
                    ans.arr[i][j] = (ans.arr[i][j] + (arr[i][k] * o.arr[k][j]) % mod) % mod;
                }
            }
        }
        return ans;
    }

    Matrix fpow(ll e)
    {
        Matrix temp = *this;
        Matrix ans;

        for (int i = 1; i <= n; i++)
            ans.arr[i][i] = 1;

        for (; e; e >>= 1, temp = temp * temp)
        {
            if (e & 1)
                ans = ans * temp;
        }

        return ans;
    }
};

struct Query
{
    ll k;
    ll id;

    bool operator<(const Query &o) const
    {
        return k < o.k;
    }
} ques[120];

ll ANS[120];

int main()
{
    n = 3;
    int T;
    scanf("%d", &T);

    for (int i = 0; i < T; i++)
    {
        scanf("%lld", &ques[i].k);
        ques[i].id = i;
    }

    sort(ques, ques + T);
    Matrix curr;
    ll currE = 0;
    curr.arr[1][1] = 1;
    curr.arr[2][2] = 1;
    curr.arr[3][3] = 1;
    for (int i = 0; i < T; i++)
    {
        ll k = ques[i].k;
        if (k <= 3)
        {
            ANS[ques[i].id] = 1;
            continue;
        }

        Matrix m;
        m.arr[1][1] = 1;
        m.arr[1][3] = 1;
        m.arr[2][1] = 1;
        m.arr[3][2] = 1;

        curr = curr * m.fpow(k - 3 - currE);
        currE = k - 3;

        ANS[ques[i].id] = (curr.arr[1][1] + (curr.arr[1][2] + curr.arr[1][3]) % mod) % mod;
    }

    for (int i = 0; i < T; i++)
    {
        printf("%lld\n", ANS[i]);
    }

    return 0;
}

高斯消元与线性方程、行列式与矩阵的逆

行列式与矩阵的逆

首先我们定义矩阵伴随因子也称余子式我为：

$$C[i,j]=(-1{)}^{i+j}\det (M[i,j])$$

其中$M[i,j]$是原矩阵$A$去掉$i$行和$j$列后剩下的矩阵。

则行列式的定义（按照第一行展开）为：

$$\det (A)=\sum _{j=1}^{n}A[1,j]C[1,j]$$

还可以按照某一行某一列展开均可。

我们定义矩阵$A$的逆矩阵为：

$$A{A}^{-1}=I$$

并且

$$A^{-1}[i,j]=\frac{C[j,i]}{\det A}$$

一个方阵有逆矩阵当且仅当这个矩阵$\det A\ne 0$。

P4783

高斯消元与线性方程组

每一个线性方程组都可以写成是一个矩阵方程$A\vec{x}=\vec{b}$，其中$A$称为系数矩阵，$\vec{x}$为解向量，$\vec{b}$为常数向量。

高斯消元法理论的核心主要如下：

• 两方程互换，解不变；
• 方程乘以非零数$k$ ，解不变；
• 方程乘以数$k$加上另一方程，解不变。

高斯消元法在将增广矩阵化为最简形后对于自由未知量的赋值，需要掌握线性相关知识，且赋值存在人工经验的因素，使得在学习过程中有一定的困难，将高斯消元法划分为五步骤，从而提出五步骤法，内容如下：

1. 增广矩阵行初等行变换为行最简形；
2. 还原线性方程组；
3. 求解第一个变量；
4. 补充自由未知量；
5. 列表示方程组通解。

P3389 弱-高斯消元
P2455 强-高斯消元

#include <bits/stdc++.h>

using namespace std;

#define FR freopen("in.txt", "r", stdin)
#define FW freopen("out.txt", "w", stdout)

typedef long long ll;

double matrix[105][105];
double esp = 1e-3;

int main()
{
    int n;
    scanf("%d", &n);

    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j <= n; j++)
        {
            scanf("%lf", &matrix[i][j]);
        }
    }

    for (int r = 0; r < n; r++)
    {
        int domain = r;

        // 选择系数最大的那行作为主元，为了减少精度损失
        for (int j = r; j < n; j++)
        {
            if (abs(matrix[j][r]) > abs(matrix[domain][r]))
                domain = j;
        }

        // 没有主元，跳过
        if (abs(matrix[domain][r]) < esp)
        {
            continue;
        }

        // 交换两行

        for (int j = 0; j <= n; j++)
        {
            swap(matrix[domain][j], matrix[r][j]);
        }

        // 消元

        for (int l = 0; l < n; l++)
        {
            if (l != r)
                for (int j = n; j >= r; j--)
                {
                    matrix[l][j] -= matrix[r][j] / matrix[r][r] * matrix[l][r];
                }
        }
    }

    // 检查解的存在性

    for (int r = 0; r < n; r++)
    {
        // 寻找这一行的主元
        int ptr = r;
        while (abs(matrix[r][ptr]) < esp && ptr <= n)
        {
            ptr++;
        }
        // 如果主元在常数列上，那么无解
        if (ptr == n)
        {
            printf("-1");
            return 0;
        }
    }

    // 如果存在缺少主元的情况，那么有无穷多解
    for (int r = 0; r < n; r++)
    {
        if (abs(matrix[r][r]) < esp)
        {
            printf("0");
            return 0;
        }
    }

    // 否则有唯一解
    for (int r = 0; r < n; r++)
    {
        // 除以主元
        printf("x%d=%.2lf\n", r + 1, matrix[r][n] / matrix[r][r]);
    }
    return 0;
}

通过求解矩阵的逆可以进行解线性方程组$\vec{x}=A^{-1}\vec{b}$。

具体应用

• 计算几何：通过向量、叉乘等运算去计算几何之间的关系
• 概率论：马尔科夫链、计数、矩阵快速幂优化等
• 多项式：多项式算法的基础