矩阵的算法
关于矩阵的算法,首先要介绍高斯消元。
在算法导论第三版中译本第484页上,写了高斯消元算法。但是问题在于,它的写法基于LUP分解,即下三角阵L、上三角阵U和置换阵P的分解方式解释高斯消元,于是我们看得一头雾水。这是显然的结果。
一种较好的办法,是用初等行变换来思考问题,放弃书上LUP分解的思路。下面的代码是解方程Ax1=x0,解完之后x0自动变为x1,即原结果不保存。因为,向量x也随着初等行变换一起操作。
先看代码:
void Gauss()//equ是一列的元素,var是一行的元素数
{
long long k;//k指向当前行或列
for(k=0;kequkvar;k++)
{
long long maxr=k;
long long i;//第一维,纵向
for(i=k+1;iequ;i++)
{
if(abs(A[i][k])abs(A[maxr][k]))
{
maxr=i;
}
}
if(k!=maxr)//换成绝对值最大的元素
{
long long j;//第二维,横向
for(j=k;jvar;j++)
{
swap(A[k][j],A[maxr][j]);//swap是简写,可以手写实现
}
swap(x[k],x[maxr]);
}
//以上为初等行变换:交换两行
x[k]/=A[k][k];//这里如果A[k][k]是0,代表方阵不满秩,可以写一个报错
long long j;
for(j=k+1;jvar;j++)
{
A[k][j]/=A[k][k];//这一行同时除以A[k][k]
}
A[k][k]=1;
//以上为初等行变换:某行乘非零数
for(i=0;iequ;i++)//从上到下的每一行,减去若干倍的当前行
{
if(i!=k)//只要不是这行
{
x[i]-=x[k]*A[i][k];//就把这项减掉
for(j=k+1;jvar;j++)//并把当前行
{
A[i][j]-=A[k][j]*A[i][k];//每个元素都减掉该项
}
A[i][k]=0;//这位置清零了
}
}
//以上为初等行变换:某行加上若干倍另一行
}
}
执行完这份代码之后,A变为单位阵,向量x也完成变换。
高斯消元也可以在取模意义下进行,只需要将除法替换为求逆,并且每次运算后都取模就行了。
我们可以看到,所谓的高斯消元算法,无非是在对元素循环的每一次中,做了三步:
A交换两行,找到系数绝对值最大的元素。检查对角线当前元素非0。
B将该行乘非0数,使得对角线当前元素系数变为1。
C所有行(除了当前行)减去若干倍当前行,消去此元。
每一步都是一个初等行变换。
算法导论上放弃了行变换的一致性,改用LUP分解:
交换两行的操作,对应置换阵P。
该行乘非零数这步和第三步合并了,变为:
For i=k+1 to n
A[i][k]/=A[k][k]
For j=k+1 to n
A[i][j]-=A[i][k]*A[k][j]
即,只对对角线当前元素的右下角子矩阵进行操作:对角线元素下方的部分进行乘非0数操作,对角线右下方的部分做减去若干倍当前行的操作。
这样一来,就将方阵A变为了PA=LU,L是对角线均为1的下三角阵,U是上三角阵。算法结束后,A变为L和U的拼接,对角线上为U。
这就导致算法导论上的算法华而不实。
以上利用初等行变换,已经解出了线性方程Ax=b。下一个问题,如何对方阵求逆?
上面的算法是,如果对方阵A进行初等行变换,变为了单位阵I,那么对于向量b执行完全相同的初等行变换,就会变为需要求解的向量x。
一个显然的事实是,如果对方阵A进行初等行变换,变为了单位阵I,那么对于单位阵I执行完全相同的初等行变换,就会变为方阵A的逆。
因此,只需要把上文的向量x替换成单位阵I就行了。以下假定算法执行前,I已经被初始化为单位阵。
void Inverse ()
{
long long k;
for(k=0;kequkvar;k++)
{
long long maxr=k;
long long i;
for(i=k+1;iequ;i++)
{
if(abs(A[i][k])abs(A[maxr][k]))
{
maxr=i;
}
}
if(k!=maxr)
{
long long j;
for(j=k;jvar;j++)
{
swap(A[k][j],A[maxr][j]);
}
long long t;
for(t=0;tvar;t++)//交换两行
{
swap(I[k][t],I[maxr][t]);
}
}
long long t;
for(t=0;tvar;t++)//乘非0数
{
I[k][t]/=A[k][k];
}
long long j;
for(j=k+1;jvar;j++)
{
A[k][j]/=A[k][k];
}
A[k][k]=1;
for(i=0;iequ;i++)
{
if(i!=k)
{
for(t=0;tvar;t++)//加某行若干倍
{
I[i][t]-=I[k][t]*A[i][k];
}
for(j=k+1;jvar;j++)
{
A[i][j]-=A[k][j]*A[i][k];
}
A[i][k]=0;
}
}
}
}
为了方便理解,这里用一个新游标t来跑单位阵I的一行。事实上完全可以用j代替。
那么行列式的办法也就清楚了。上述操作中,只有某行乘非0数会改变行列式的值,交换两行行列式取负,那么在这里修改行列式值就行了。如果不满秩,行列式为0。
void Det()
{
long long det=1;
long long k;
for(k=0;kequkvar;k++)
{
long long maxr=k;
long long i;
for(i=k+1;iequ;i++)
{
if(abs(A[i][k])abs(A[maxr][k]))
{
maxr=i;
}
}
if(k!=maxr)
{
long long j;
for(j=k;jvar;j++)
{
swap(A[k][j],A[maxr][j]);
}
det*=-1;
}
if(A[k][k]==0)
{
return 0;
}
det*=A[k][k];
long long j;
for(j=k+1;jvar;j++)
{
A[k][j]/=A[k][k];
}
A[k][k]=1;
for(i=0;iequ;i++)
{
if(i!=k)
{
for(j=k+1;jvar;j++)
{
A[i][j]-=A[k][j]*A[i][k];
}
A[i][k]=0;
}
}
}
return det;
}
另一个与矩阵相关的算法是最小二乘法。
对于一个给定的关于自变量x的矩阵函数A(x),希望找到一个解向量c,使得对于给定的结果向量y,||A(x)c-y||最小。符号||x||,表示向量x乘向量x的转置再开二次根号,即它与它自己的内积的平方根,即向量模长。
因为矩阵函数A(x)是已知的,是设计拟合函数的时候已经设计好的。比如,A(x)是幂函数,那么A(x)乘解向量c得到的A(x)c就是多项式。那么,由于样本x已知,直接带入A(x),就能计算出已知的矩阵A。当然,对于样本x,给定的结果y也是已知的。上文当中A和y都已知,只有待求的c不知道。
显然。如果Ac=y有解,自然而然,差是零向量,那么向量模长当然最小,用上文的高斯消元,问题就解决了。但是,我们不知道它有没有解,有没有解要高斯消元跑完才知道,并且高斯消元解决不了无解的情况。
这时,就需要最小二乘法。对||A(x)c-y||模长平方,展开,求关于c的微分。令微分为0,得到方程:
AT(Ac-y)=0
A的转置乘(Ac-y)为0,即ATAc=ATy。(T写在右上角,表示转置)
只要A是列满秩,则ATA是可逆的矩阵。当然,ATA是对称阵,A列满秩时ATA是正定阵。我们可以把A想象为一个瘦高的矩阵,即样本(x,y)足够多,使得拟合容易进行。
A的列满秩,也可以理解为A的各个列之间不相关,即要求的系数c之间不相关,最终拟合的函数各个项之间不相关。
于是,只需要利用上述算法解方程ATAc=ATy就行了,在ATA可逆的时候一定有解。
稀疏矩阵算法是什么
稀疏矩阵算法是以稀疏矩阵作为核心数据结构的算法。
稀疏矩阵算法的最大特点是通过只存储和处理非零元素从而大幅度降低存储空间需求以及计算复杂度,代价则是必须使用专门的稀疏矩阵压缩存储数据结构。稀疏矩阵算法是典型的不规则算法,计算访存比很低,并且计算过程中的访存轨迹与稀疏矩阵的稀疏结构相关。
稀疏矩阵算法是自然科学和社会科学中许多领域进行数值模拟计算时的关键技术和性能瓶颈,为了提高稀疏矩阵算法的计算性能,需要提高相应算法在特定平台上的计算效率。
矩阵的逆矩阵怎么求
初等行变换不影响线性方程组的解,也可用于高斯消元法,用于逐渐将系数矩阵化为标准形。初等行变换不改变矩阵的核(故不改变解集),但改变了矩阵的像。反过来,初等列变换没有改变像却改变了核。
矩阵的逆矩阵怎么求
运用初等行变换法。将一n阶可逆矩阵A和n阶单位矩阵I写成一个nX2n的矩阵B=(A,I])对B施行初等行变换,即对A与I进行完全相同的若干初等行变换,目标是把A化为单位矩阵。当A化为单位矩阵I的同时,B的右一半矩阵同时化为了A的逆矩阵。
逆矩阵的性质
可逆矩阵一定是方阵。
如果矩阵A是可逆的,其逆矩阵是唯一的。
A的逆矩阵的逆矩阵还是A。记作(A-1)-1=A。
可逆矩阵A的转置矩阵AT也可逆,并且(AT)-1=(A-1)T (转置的逆等于逆的转置)。
若矩阵A可逆,则矩阵A满足消去律。即AB=O(或BA=O),则B=O,AB=AC(或BA=CA),则B=C。
两个可逆矩阵的乘积依然可逆。
矩阵可逆当且仅当它是满秩矩阵。
|矩阵的算法
矩阵的算法 矩阵的逆矩阵怎么求 稀疏矩阵算法