算法分析与设计 - 作业9

目录

• 问题一
• 问题二

问题一

调研学习并给出矩阵的LU分解方法。

\(LU\) 分解是“矩阵因式分解”的一种，旨在将某个矩阵 \(A\) 表示为两个或多个矩阵的乘积。具体地，\(LU\) 分解可将矩阵表示为 \(A=LU\) 的形式，其中矩阵 \(L\) 代表Lower Triangular（下三角矩阵），矩阵 \(U\) 代表Upper Triangular（上三角矩阵）。如：

\[A = \begin{pmatrix} 1& 1 &-1\\ 1& 2 &1\\ -1& 1 &6 \end{pmatrix} = \begin{pmatrix} 1& 0 &0\\ 1& 1 &0\\ -1& 2 &1 \end{pmatrix} \begin{pmatrix} 1& 1 &-1\\ 0& 1 &2\\ 0& 0 &1 \end{pmatrix}\]

能够 \(LU\) 分解的矩阵需要满足以下三个条件：

• 矩阵是方阵；
• 矩阵是可逆的，也就是该矩阵是满秩矩阵，每一行都是独立向量；
• 消元过程中没有 0 主元出现，也就是消元过程中不能出现行交换的初等变换。

进行 LU 分解可考虑按照矩阵乘法的原则进行递推。即有下式：

\[\begin{cases} &u_{1i} =a_{1i} &(i=1,2,\cdots ,n)\\ &l_{i1} =a_{i1}/u_{11} &(i=2,3,\cdots ,n)\\ &u_{ri} =a_{ri}-\sum_{k=1}^{r-1}l_{rk}u_{ki} &(i=r,r+1,\cdots ,n)\\ &l_{ir} =\dfrac{(a_{ir}-\sum_{k=1}^{r-1}l_{ik}u_{kr})}{u_{rr}} &((i=r+1,\cdots ,n) \land (r\ne n)) \end{cases}\]

代码：

#include <iostream>
// 参数：一个order阶矩阵，和矩阵的阶数
void lowerUpperFactor(double *matrix, int order) {
    printf("--------原矩阵：--------\n");
    printMatrix(matrix, order,order);
    
    // 结果变量  L矩阵和U矩阵都是order阶矩阵
    double *L = new double[order*order];
    double *U = new double[order*order];
    // 初始化全为0
    for (int i = 0; i < order; i++) {
        // 初始化U下三角为0
        for (int j = 0; j < i; j++) {
            *(U + i * order + j) = 0;
        }
        //初始化L对角线为1，上三角为0
        *(L + i * order + i) = 1;
        for (int j = i + 1; j < order; j++) {
            *(L + i * order + j) = 0;
        }
    }
    // 计算U的第一行和L的第一列
    int i = 0;
    for (i = 0; i < order; i++) {
        *(U + i) = *(matrix + i);
    }
    for (i = 1; i < order; i++) {
        *(L + i * order) = *(matrix + i * order) / *U;
    }
    // 计算其余行列
    int temp;
    for (int i = 1; i < order; i++) {
        // 计算矩阵U
        for (int j = i; j < order; j++) {
            temp = 0;
            for (int k = 0; k < i; k++) {
                temp+= (*(U + k * order + j) * (*(L + i * order + k)));
            }
            *(U + i * order + j) = *(matrix + i * order + j) - temp;
        }
        // 计算矩阵L
        for (int j = i+1; j < order; j++) {
            temp = 0;
            for (int k = 0; k < i; k++) {
                temp += *(U + k * order + i) * (*(L + j * order + k));
            }
            *(L + j * order + i) = (*(matrix +j * order + i) - temp) / (*(U+i* order + i));
        }
    }

    printf("------矩阵U------\n");
    printMatrix(U, order,order);
    printf("------矩阵L------\n");
    printMatrix(L, order, order);

    if (L) {
        delete[] L;
    }
    if (U) {
        delete[] U;
    }
}
void printMatrix(double *matrix, int row, int column) {　　for (int i = 0; i < row; i++) {　　　　for (int j = 0; j < column; j++) {　　　　　　printf("%6.2lf    ", *(matrix + i * column + j));　　　　}　　　　printf("\n");　　}}
int main() {

    //double matrix[] = { 1,2,3,2,5,2,3,1,5 };
    //int order = 3;
    double matrix1[] = { 1,1,-1,2,1,2,0,2,-1,-1,2,0,0,0,-1,1 };
    int order1 = 4;

    lowerUpperFactor(matrix1, order1);
}

问题二

给出方案计算可逆矩阵的逆。

一种使用高斯-约旦消元的做法。

记 \([AB]\) 表示将一个 \(n\times m_b\) 的矩阵 \(B\) 拼接一个 \(n\times m_a\) 的矩阵 \(A\) 之后形成的矩阵，即有：

\[[AB]_{i, j} = \begin{cases} A_{i, j} &(j\le m_a)\\ B_{i, j - m_a} &(m_a <j \le m_a+m_b) \end{cases}\]

记 \(I\) 为 \(n\times n\) 的单位矩阵，则有：

\[A^{-1}\times [AI] = [IA^{-1}] \]

则若可以通过对矩阵 \([AI]\) 进行若干操作使得矩阵的左半部分（即拼接后左半边的 \(A\)）转化为单位矩阵 \(I\)，则我们进行的操作等价于令 \(A^{-1}\) 乘上了矩阵 \([AI]\)，则得到的矩阵的右半部分即为所求的逆矩阵 \(A^{-1}\)。

使用高斯-约旦法将 \([AI]\) 的左半边消成对角矩阵，再令每行除以系数即可。

下列代码在模意义下实现了矩阵求逆，使用逆元替换了消元时的除法。

总复杂度 \(O(n^3)\) 级别。

代码：

//By:Luckyblock
/*
*/
#include <cstdio>
#include <cctype>
#include <algorithm>
#define LL long long
const int kN = 410;
const LL p = 1e9 + 7;
//=============================================================
int n;
LL a[kN][kN << 1];
//=============================================================
inline int read() {
  int f = 1, w = 0; char ch = getchar();
  for (; !isdigit(ch); ch = getchar()) if (ch == '-') f = - 1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + ch - '0';
  return f * w;
}
LL qpow(LL x_, LL y_) {
  LL ret = 1;
  while (y_) {
    if (y_ & 1) ret = ret * x_ % p;
    x_ = x_ * x_ % p, y_ >>= 1ll;
  }
  return ret;
}
bool Gauss_jordan() {
  for (int i = 1; i <= n; ++ i) {
    int pos = i;
    for (int j = i + 1; j <= n; ++ j) {
      if (a[j][i] > a[pos][i]) pos = j;
    }
    if (pos != i) std::swap(a[i], a[pos]);
    if (!a[i][i]) return false;

    LL inv = qpow(a[i][i], p - 2);
    for (int j = 1; j <= n; ++ j) {
      if (j == i) continue;
      LL delta = a[j][i] * inv % p;
      for (int k = 1; k <= 2 * n; ++ k) {
        a[j][k] = ((a[j][k] - delta * a[i][k]) % p + p) % p;
      }
    }
    for (int j = 1; j <= 2 * n; ++ j) a[i][j] = a[i][j] * inv % p;
  }
  return true;
}
//=============================================================
int main() {
  // freopen("1.txt", "r", stdin);
  n = read();
  for (int i = 1; i <= n; ++ i) {
    for (int j = 1; j <= n; ++ j) {
      a[i][j] = read(), a[i][i + n] = 1;
    }
  }
  if (!Gauss_jordan()) {
    printf("No Solution\n");
    return 0;
  }
  for (int i = 1; i <= n; ++ i) {
    for (int j = n + 1; j <= 2 * n; ++ j) {
      printf("%lld ", a[i][j]);
    }
    printf("\n");
  }
  return 0;
}