算法导论28（矩阵运算）

最新推荐文章于 2022-07-26 19:02:24 发布

原创最新推荐文章于 2022-07-26 19:02:24 发布 · 1.5k 阅读

0 ·

CC 4.0 BY-SA版权

算法导论专栏收录该内容

22 篇文章

订阅专栏

本文深入探讨了线性方程组的求解方法，包括LU分解、矩阵求逆及对称正定矩阵的应用。通过实例演示了如何使用LU分解法求解线性方程组，并介绍了最小二乘逼近的概念。此外，还详细解释了矩阵求逆的过程，以及如何通过矩阵求逆解决特定问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

28.1 求解线性方程组
$Ax = b,PA = LU \Rightarrow LUx = Pb \Rightarrow Ly = Pb,Ux = y$

28.2 矩阵求逆

A X = E \Rightarrow A x i = e i A x i = e i, P A = L U \Rightarrow L y i = P e i, U x i = y i

$\begin{array}{l} AX = E \Rightarrow A{x_i} = {e_i}\\ A{x_i} = {e_i},PA = LU \Rightarrow L{y_i} = P{e_i},U{x_i} = {y_i} \end{array}$

28.3 对称正定矩阵和最小二乘逼近

推论28.6：一个对称正定矩阵的 $LU$ 分解永远不会出现除数为0的情形。

最小二乘逼近
${A^T}Ac = {A^T}y \Rightarrow c = {({A^T}A)^{ - 1}}{A^T}y \Rightarrow {A^ + } = {({A^T}A)^{ - 1}}{A^T},c = {A^ + }y$

vector<double> LUPSolve(vector<vector<double> > L,vector<vector<double> > U,vector<double> pi,vector<double> b)
{
    int n=L.size();
    vector<double> x(n,0),y(n,0);
    //Ly=Pb，正向替换求y
    for(int i=0;i<n;++i)
    {
        y[i]=b[pi[i]];
        for(int j=0;j<i;++j)y[i]-=L[i][j]*y[j];
    }
    //Ux=y，反向替换求x
    for(int i=n-1;i>=0;--i)
    {
        x[i]=y[i];
        for(int j=i+1;j<n;++j)x[i]-=U[i][j]*x[j];
        x[i]/=U[i][i];
    }
    return x;
}

//LU分解
void LUDecomposition(vector<vector<double> > A,vector<vector<double> > &L,vector<vector<double> > &U)
{
    int n=A.size();
    for(int i=0;i<n;++i)L[i][i]=1;
    for(int i=0;i<n;++i)
    {
        for(int j=i+1;j<n;++j)L[i][j]=0;
    }
    for(int i=0;i<n;++i)
    {
        for(int j=0;j<i;++j)U[i][j]=0;
    }
    for(int k=0;k<n;++k)
    {
        U[k][k]=A[k][k];
        for(int i=k+1;i<n;++i)
        {
            L[i][k]=A[i][k]/U[k][k];
            U[k][i]=A[k][i];
        }
        for(int i=k+1;i<n;++i)
        {
            for(int j=k+1;j<n;++j)A[i][j]-=L[i][k]*U[k][j];
        }
    }
}

//LUP分解
bool LUPDecomposition(vector<vector<double> > &A,vector<double> &pi)
{
    int n=A.size();
    for(int i=0;i<n;++i)pi[i]=i;
    for(int k=0;k<n;++k)
    {
        int p=0,k0=-1;
        for(int i=k;i<n;++i)
        {
            if(abs(A[i][k])>p)
            {
                p=abs(A[i][k]);
                k0=i;
            }
        }
        if(p==0)return false;
        swap(pi[k],pi[k0]);
        for(int i=0;i<n;++i)swap(A[k][i],A[k0][i]);
        for(int i=k+1;i<n;++i)
        {
            A[i][k]/=A[k][k];
            for(int j=k+1;j<n;++j)A[i][j]-=A[i][k]*A[k][j];
        }
    }
    return true;
}