矩阵乘法是线性代数中最常见的运算之一,它在数值计算中有广泛的应用。若A和B是2个nn的矩阵,则它们的乘积C=AB同样是一个nn的矩阵。A和B的乘积矩阵C中的元素C[i,j]定义为:
若依此定义来计算A和B的乘积矩阵C,则每计算C的一个元素C[i,j],需要做n个乘法和n-1次加法。因此,求出矩阵C的n∧2个元素所需的计算时间为0(n∧3)。
20世纪60年代末,Strassen采用了类似于在大整数乘法中用过的分治技术,将计算2个n阶矩阵乘积所需的计算时间改进到O(n∧log7)=O(n∧2.18),其基本思想还是使用分治法。
首先,我们还是需要假设n是2的幂。将矩阵A,B和C中每一矩阵都分块成为4个大小相等的子矩阵,每个子矩阵都是(n/2)*(n/2)的方阵。由此可将方程C=AB重写为:
C∨11 C∨12 A∨11 A∨12 B∨11 B∨12
C∨21 C∨22 = A∨21 A∨22 = B∨21 B∨22
由此可得:
C∨11=A∨11B∨11+A∨12B∨21 (2)
C∨12=A∨11B∨12+A∨12B∨22 (3)
C∨21=A∨21B∨11+A∨22B∨21 (4)
C∨22=A∨21B∨12+A∨22B∨22 (5)
如果n=2,则2个2阶方阵的乘积可以直接用(2)-(3)式计算出来,共需8次乘法和4次加法。当子矩阵的阶大于2时,为求2个子矩阵的积,可以继续将子矩阵分块,直到子矩阵的阶降为2。这样,就产生了一个分治降阶的递归算法。依此算法,计算2个n阶方阵的乘积转化为计算8个n/2阶方阵的乘积和4个n/2阶方阵的加法。2个n/2n/2矩阵的加法显然可以在c*n2/4时间内完成,这里c是一个常数。因此,上述分治法的计算时间耗费T(n)应该满足:
这个递归方程的解仍然是T(n)=O(n3)。因此,该方法并不比用原始定义直接计算更有效。究其原因,乃是由于式(2)-(5)并没有减少矩阵的乘法次数。而矩阵乘法耗费的时间要比矩阵加减法耗费的时间多得多。要想改进矩阵乘法的计算时间复杂性,必须减少子矩阵乘法运算的次数。按照上述分治法的思想可以看出,要想减少乘法运算次数,关键在于计算2个2阶方阵的乘积时,能否用少于8次的乘法运算。Strassen提出了一种新的算法来计算2个2阶方阵的乘积。他的算法只用了7次乘法运算,但增加了加、减法的运算次数。这7次乘法是:
M∨1=A∨11(B∨12-B∨22)
M∨2=(A∨11+A∨12)B∨22
M∨3=(A∨21+A∨22)B∨11
M∨4=A∨22(B∨21-B∨11)
M∨5=(A∨11+A∨22)(B∨11+B∨22)
M∨6=(A∨12-A∨22)(B∨21+B∨22)
M∨7=(A∨11-A∨21)(B∨11+B∨12)
做了这7次乘法后,再做若干次加、减法就可以得到:
C∨11=M∨5+M∨4-M∨2+M∨6
C∨12=M∨1+M∨2
C∨21=M∨3+M∨4
C∨22=M∨5+M∨1-M∨3-M∨7
以上计算的正确性很容易验证。例如:
C∨22=M∨5+M∨1-M∨3-M∨7
=(A∨11+A∨22)(B∨11+B∨22)+A∨11(B∨12-B∨22)-(A∨21+A∨22)B∨11-(A∨11-A∨21)(B∨11+B∨12)
=A∨11B∨11+A∨11B∨22+A∨22B∨11+A∨22B∨22+A∨11B∨12
-A∨11B∨22-A∨21B∨11-A∨22B∨11-A∨11B∨11-A∨11B∨12+A∨21B∨11+A∨21B∨12
=A∨21B∨12+A∨22B∨22
Strassen矩阵乘积分治算法中,用了7次对于n/2阶矩阵乘积的递归调用和18次n/2阶矩阵的加减运算。由此可知,该算法的所需的计算时间T(n)满足如下的递归方程:
按照解递归方程的套用公式法,其解为T(n)=O(nlog7)≈O(n2.81)。由此可见,Strassen矩阵乘法的计算时间复杂性比普通矩阵乘法有阶的改进。
有人曾列举了计算2个2阶矩阵乘法的36种不同方法。但所有的方法都要做7次乘法。除非能找到一种计算2阶方阵乘积的算法,使乘法的计算次数少于7次,按上述思路才有可能进一步改进矩阵乘积的计算时间的上界。但是Hopcroft和Kerr(197l)已经证明,计算2个22矩阵的乘积,7次乘法是必要的。因此,要想进一步改进矩阵乘法的时间复杂性,就不能再寄希望于计算22矩阵的乘法次数的减少。或许应当研究33或55矩阵的更好算法。在Strassen之后又有许多算法改进了矩阵乘法的计算时间复杂性。目前最好的计算时间上界是O(n2.367)。而目前所知道的矩阵乘法的最好下界仍是它的平凡下界Ω(n2)。因此到目前为止还无法确切知道矩阵乘法的时间复杂性。关于这一研究课题还有许多工作可做。
Demo:
package Recursive;
/**
* Strassen 矩形算法
*/
import java.io.*;
import java.util.*;
class matrix
{
public int[][] m = new int[32][32];
}
public class Strassen
{
public int judgment(int n)
{
int flag = 0,temp=n;
while(temp%2==0)
{
if(temp%2==0)
temp/=2;
else flag=1;
}
if(temp==1)
flag=0;
return flag;
}
public void Divide(matrix d,matrix d11,matrix d12,matrix d21,matrix d22,int n)/**//*分解矩阵方法*/
{
int i,j;
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
{
d11.m[i][j]=d.m[i][j];
d12.m[i][j]=d.m[i][j+n];
d21.m[i][j]=d.m[i+n][j];
d22.m[i][j]=d.m[i+n][j+n];
}
}
public matrix Merge(matrix a11,matrix a12,matrix a21,matrix a22,int n)/**//*合并矩阵方法*/
{
int i,j;
matrix a = new matrix();
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
{
a.m[i][j]=a11.m[i][j];
a.m[i][j+n]=a12.m[i][j];
a.m[i+n][j]=a21.m[i][j];
a.m[i+n][j+n]=a22.m[i][j];
}
return a;
}
public matrix AdhocMatrixMultiply(matrix x,matrix y) /**//*阶数为2的矩阵乘法方法*/
{
int m1,m2,m3,m4,m5,m6,m7;
matrix z = new matrix();
m1=(y.m[1][2] - y.m[2][2]) * x.m[1][1];
m2=y.m[2][2] * (x.m[1][1] + x.m[1][2]);
m3=(x.m[2][1] + x.m[2][2]) * y.m[1][1];
m4=x.m[2][2] * (y.m[2][1] - y.m[1][1]);
m5=(x.m[1][1] + x.m[2][2]) * (y.m[1][1]+y.m[2][2]);
m6=(x.m[1][2] - x.m[2][2]) * (y.m[2][1]+y.m[2][2]);
m7=(x.m[1][1] - x.m[2][1]) * (y.m[1][1]+y.m[1][2]);
z.m[1][1] = m5 + m4 - m2 + m6;
z.m[1][2] = m1 + m2;
z.m[2][1] = m3 + m4;
z.m[2][2] = m5 + m1 - m3 - m7;
return z;
}
public matrix MatrixPlus(matrix f,matrix g,int n) /**//*矩阵加法方法*/
{
int i,j;
matrix h = new matrix();
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
h.m[i][j]=f.m[i][j]+g.m[i][j];
return h;
}
public matrix MatrixMinus(matrix f,matrix g,int n) /**//*矩阵减法方法*/
{
int i,j;
matrix h = new matrix();
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
h.m[i][j]=f.m[i][j]-g.m[i][j];
return h;
}
public matrix MatrixMultiply(matrix a,matrix b,int n) /**//*矩阵乘法方法*/
{
int k;
matrix a11,a12,a21,a22;
a11 = new matrix();
a12 = new matrix();
a21 = new matrix();
a22 = new matrix();
matrix b11,b12,b21,b22;
b11 = new matrix();
b12 = new matrix();
b21 = new matrix();
b22 = new matrix();
matrix c11,c12,c21,c22,c;
c11 = new matrix();
c12 = new matrix();
c21 = new matrix();
c22 = new matrix();
c = new matrix();
matrix m1,m2,m3,m4,m5,m6,m7;
k=n;
if(k==2)
{
c=AdhocMatrixMultiply(a,b);
return c;
}
else
{
k=n/2;
Divide(a,a11,a12,a21,a22,k); //拆分A、B、C矩阵
Divide(b,b11,b12,b21,b22,k);
Divide(c,c11,c12,c21,c22,k);
m1=MatrixMultiply(MatrixMinus(b12,b22,n/2),a11,k);
m2=MatrixMultiply(b22,MatrixPlus(a11,a12,k),k);
m3=MatrixMultiply(MatrixPlus(a21,a22,k),b11,k);
m4=MatrixMultiply(a22,MatrixMinus(b21,b11,k),k);
m5=MatrixMultiply(MatrixPlus(a11,a22,k),MatrixPlus(b11,b22,k),k);
m6=MatrixMultiply(MatrixMinus(a12,a22,k),MatrixPlus(b21,b22,k),k);
m7=MatrixMultiply(MatrixMinus(a11,a21,k),MatrixPlus(b11,b12,k),k);
c11=MatrixPlus(MatrixMinus(MatrixPlus(m5,m4,k),m2,k),m6,k);
c12=MatrixPlus(m1,m2,k);
c21=MatrixPlus(m3,m4,k);
c22=MatrixMinus(MatrixMinus(MatrixPlus(m5,m1,k),m3,k),m7,k);
c=Merge(c11,c12,c21,c22,k); //合并C矩阵
return c;
}
}
public static void main(String[] args)throws IOException
{
Strassen instance = new Strassen();
int i,j,num;
matrix A,B,C;
A = new matrix();
B = new matrix();
C = new matrix();
Scanner in = new Scanner(System.in);
System.out.print("输入矩阵的阶数: ");
num = in.nextInt();
if(instance.judgment(num)==0)
{
System.out.println("输入矩阵A:");
for(i=1;i<=num;i++)
for(j=1;j<=num;j++)
A.m[i][j] = in.nextInt();
System.out.println("输入矩阵B:");
for(i=1;i<=num;i++)
for(j=1;j<=num;j++)
B.m[i][j] = in.nextInt();
if(num==1)
C.m[1][1]=A.m[1][1]*B.m[1][1]; //矩阵阶数为1时的特殊处理
else
C=instance.MatrixMultiply(A,B,num);
System.out.println("矩阵C为:");
for(i=1;i<=num;i++)
{
for(j=1;j<=num;j++)
System.out.print(C.m[i][j] + " ");
System.out.println();
}
}
else
System.out.println("输入的阶数不是2的N次方");
}
}
显示结果为
输入矩阵的阶数: 2
输入矩阵A:
12
12
12
12
输入矩阵B:
12
12
12
12
矩阵C为:
288 288
288 288