算法导论28（矩阵运算）

28.1 求解线性方程组

Ax=b,PA=LU⇒LUx=Pb⇒Ly=Pb,Ux=yAx = b,PA = LU \Rightarrow LUx = Pb \Rightarrow Ly = Pb,Ux = y

28.2 矩阵求逆

AX=E⇒Axi=eiAxi=ei,PA=LU⇒Lyi=Pei,Uxi=yi\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：一个对称正定矩阵的LULU分解永远不会出现除数为0的情形。

最小二乘逼近

ATAc=ATy⇒c=(ATA)−1ATy⇒A+=(ATA)−1AT,c=A+y{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;
}