title: 稀疏矩阵乘法
 date: 2020-11-09 19:31:44
 tags: 稀疏矩阵运算
 categories: 数据结构
 
 

 在本算法中，两个稀疏矩阵的特性都有用到
 
 
规定
 
规定以下变量名称，本文讲述 矩阵A × 矩阵B = 矩阵C 的运算过程
 
需要用到的存储结构有：
 
矩阵A，矩阵 B 的原始二维数组(2个)
矩阵A，矩阵B 的三元组数组(2个)
存储 矩阵A，矩阵B 每行有多少个非零元的数组(2个，分别存A、B矩阵)
存储 矩阵B每行首个非零元在三元组数组中的位置的数组(1个)
需要开一个中间数组用于存矩阵C的每一行 （这个中间数组的长度等于B的列数，它只负责在一次循环中记录矩阵C的原始的一行，注意不是稀疏矩阵表示的一样，而是原始矩阵C的一行）
矩阵C 的三元组数组，将稀疏矩阵还原以后的 C的二维数组
 
综上，本算法需要 5 种，9 个数组
 
思路
 
首先，将上述1 - 4 中的 7 个数组均求出来
 
一切准备妥当后，且听我慢慢道来
 
遍历A的三元组数组， 将位于同一行的每个非零元分别与 该非零元的列数 对应到B中的行数，对B中行数是 非零元列数的每个非零元依次求乘法运算并存到中间数组的对应下标位置处，当遍历完A的一行的所有非零元之后，也就求出了C的一行(存在中间数组中，这个中间数组存的不是稀疏的，而是原始矩阵C的一行)，然后再对C进行压缩处理，因为在这个过程中，可能原来不是非零元，在计算完之后成了非零元，压缩完之后存到C的三元组数组中
 
代码
 
此代码仅仅是将每次内层for循环得到的temp数组(即C矩阵的每一行)打印了出来，还没有做对C矩阵的录入工作，和前面所说稍有出入，还缺个尾巴，但核心算法问题已经解决。
 
#include <stdio.h>
#include <stdlib.h>
/*
    矩阵 a    4 × 7   8个非零元
    0 8 0 0 6 0 0 
    0 1 0 3 0 0 0 
    7 0 0 2 0 4 0 
    0 0 8 0 0 0 5

    矩阵 b    7 × 3   6个非零元
    1  0  0
    0  0  2
    0  0  0
    0  9  0 
    3  0  0
    0  0  6
    0  4  0
*/

// 定义非零元的结构体(三元组的基本结构)
typedef struct member
{
    // 矩阵中非零元的 行，列，值
    int row, col, x;
    
} Member;

// 定义稀疏矩阵
typedef struct matrix
{
    // 三元组数组
    Member *arr;
    // 矩阵的行数、列数、非零元个数
    int rows, cols, counts;
    
} Matrix;

// 初始化矩阵A的三元组数组
void InitMatrixA(Matrix *);
// 初始化矩阵B的三元组数组
void InitMatrixB(Matrix *);
// 初始化一个三元组数组(在本算法中用于初始化矩阵C的压缩矩阵表示)
void InitMatrix(Matrix *);
// 以原始形式打印矩阵
void PrintMatrix(Matrix);

// 求存储了矩阵每一行的第一个非零元在三元组数组中的下标是多少的数组
int *getIndexArr(Matrix);
int *getCountArr(Matrix);
// 求每一行有多少个非零元
int *getNotZeroSuchRow(Matrix);
// 矩阵乘法
void Multiply(Matrix, Matrix, int *, int *, int *, int *);

int main(void)
{
    Matrix a, b;
    InitMatrixA(&a);
    int *arrA1 = getIndexArr(a);
    int *arrA2 = getNotZeroSuchRow(a);
    InitMatrixB(&b);
    int *arrB1 = getIndexArr(b);
    int *arrB2 = getNotZeroSuchRow(b);
    Multiply(a, b, arrA1, arrB1, arrA2, arrB2);
    return 0;
}

void InitMatrixA(Matrix *pMatrix)
{
    pMatrix->cols = 7;
    pMatrix->rows = 4;
    pMatrix->counts = 9;
    Member *p = (Member *)malloc(sizeof(Member) * pMatrix->counts + 1);
    p[1].row = 1; p[1].col = 2; p[1].x = 8;
    p[2].row = 1; p[2].col = 5; p[2].x = 6;
    p[3].row = 2; p[3].col = 2; p[3].x = 1;
    p[4].row = 2; p[4].col = 4; p[4].x = 3;
    p[5].row = 3; p[5].col = 1; p[5].x = 7;
    p[6].row = 3; p[6].col = 4; p[6].x = 2;
    p[7].row = 3; p[7].col = 6; p[7].x = 4;
    p[8].row = 4; p[8].col = 3; p[8].x = 8;
    p[9].row = 4; p[9].col = 7; p[9].x = 5;
    pMatrix->arr = p;
    printf("矩阵A:\n");
    PrintMatrix(*pMatrix);
}

void InitMatrixB(Matrix *pMatrix)
{
    pMatrix->rows = 7;
    pMatrix->cols = 3;
    pMatrix->counts = 6;
    Member *p = (Member *)malloc(sizeof(Member) * pMatrix->counts + 1);
    p[1].row = 1; p[1].col = 1; p[1].x = 1;
    p[2].row = 2; p[2].col = 3; p[2].x = 2;
    p[3].row = 4; p[3].col = 2; p[3].x = 9;
    p[4].row = 5; p[4].col = 1; p[4].x = 3;
    p[5].row = 6; p[5].col = 3; p[5].x = 6;
    p[6].row = 7; p[6].col = 2; p[6].x = 4;
    pMatrix->arr = p;
    printf("矩阵B:\n");
    PrintMatrix(*pMatrix);
}

void PrintMatrix(Matrix matrix)
{
    int index = 1;
    for (int i = 1; i <= matrix.rows; i++)
    {
        for (int j = 1; j <= matrix.cols; j++)
        {
            if (i == matrix.arr[index].row && j == matrix.arr[index].col)
                printf("%d ", matrix.arr[index++].x);
            else
                printf("0 ");
        }
        printf("\n");
    }
}

int *getIndexArr(Matrix matrix)
{
    // 要返回的目的数组
    int *p = (int *)malloc(sizeof(int) * matrix.rows);
    // 求每一行有多少个非0元的数组
    int arr[100] = {0};
    for (int i = 1; i <= matrix.counts; i++)
        arr[matrix.arr[i].row]++;
    p[1] = 1;
    for (int i = 2; i <= matrix.rows; i++)
        p[i] = p[i - 1] + arr[i - 1];
    return p;
}
int *getNotZeroSuchRow(Matrix matrix)
{
    int *p = (int *)malloc(sizeof(int) * (matrix.rows + 1));
    // 将每个数组元素初始化为0
    for (int i = 0; i < matrix.rows + 1; i++)
        p[i] = 0;
    for (int i = 1; i <= matrix.counts; i++)
        p[matrix.arr[i].row]++;
    return p;
}
void Multiply(Matrix a, Matrix b, int *arr1, int *brr1, int *arr2, int *brr2)
{
    int *temp = (int *)malloc(sizeof(int) * (b.cols + 1));
    for (int i = 0; i < b.cols + 1; i++)
        temp[i] = 0;
    int index = 1;
    printf("矩阵 C = A × B:\n");
    // arr和brr分别存储了矩阵A矩阵B每行首个非零元在三元组数组中的下标
    for (int i = 1; i <= a.rows; i++)
    {
        for (int j = 1; j <= arr2[i]; j++)
        {
            // a的三元组数组中每个元素的列数
            Member aMember = a.arr[index];
            int colOfSuchMemberA = aMember.col;
            // 去b矩阵的brr数组中找到行为上面得到的列数的第一个非零元的下标
            int indexOfNotZeroOfRowOfB = brr1[colOfSuchMemberA];
            for (int i = 1; i <= brr2[colOfSuchMemberA]; i++)
            {
                Member bMember = b.arr[indexOfNotZeroOfRowOfB + i - 1];
                temp[bMember.col] += bMember.x * aMember.x;
            }
            index++;
        }
        // 内层循环一次，代表C矩阵的一行已经被计算出来
        for (int i = 1; i <= b.cols; i++)
        {
            printf("%d ", temp[i]);
            temp[i] = 0;
        }
        printf("\n");
    }
}

行索引稀疏矩阵乘法（严蔚敏版）
 
原理
 
行索引稀疏矩阵查找某一列的元素没那么方便，所以在做矩阵乘法时（这里以M乘N=Q为例），严书的做法是：在求Q的某一行的值是，用M的一行去遍历N的每一行，对结果中同列的值进行累加，最后稀疏化存入Q的当前行中，这个过程还原成正常矩阵比较容易理解：
 求Q（2,2）的第一行时，肯定是M的第一行和N的第一列逐乘再累加，然后再M的第一行和N的第二列逐乘累加
 M（2,3）* N（3,2）：
 
 直接算的话就是:
 Q[1,1] = 1*1 + 2*0 + 3*1 = 1+0+3 = 4
 Q[1,2] = 1*0 + 2*1 + 3*1 = 0+2+3 = 5
 这回我们直接同时求两列的值，其实就是改变一下计算顺序：
 
 这个就相当于严书代码中ctemp的作用,只不过ctemp是直接累加.
 
代码
 
#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
const int MAX_E = 1e+4+5, MAX_R = 105;

//严版特色,下标从1开始
struct Triple{int i,j,e;};
struct RLSMatrix{
    Triple data[MAX_E];
    int rpos[MAX_R];//行首索引, rpos[i]指向第i行中的首元素在data中的序号, 则rpos[i+1]-1指向第i行中最后一个非0元素在N.data中的序号.
    int rows, cols, elems;
};

int ctemp[MAX_R];//Q的第i行元素结果累加器,算完一行就压缩存储进Q.data中


RLSMatrix mul(RLSMatrix M, RLSMatrix N){
    RLSMatrix Q;
    Q.rows = M.rows;Q.cols = N.cols;Q.elems = 0;
    if(M.elems * N.elems != 0){                              // Q是非零矩阵
        for(int arow = 1; arow <= M.rows; arow++){
            memset(ctemp, 0, sizeof(ctemp));              //到了下一行,清零
            Q.rpos[arow] = Q.elems+1;                  //新一行的行首索引是当前data中元素个数+1,从该处继续向data中存三元组

            int tp;      //用tp指向data中该行行末的下一个位置,方便遍历
            if(arow < M.rows) tp = M.rpos[arow+1];    //如果不是最后一行, tp指向下一行行首
            else tp = M.elems+1;   //最后一行, tp直接指向data中最后一个+1

            for(int p = M.rpos[arow];p < tp; p++){  //对当前行找到的每一个非零元
                int brow = M.data[p].j;     //找到对应元在N中的行号(M中某行第三列的元素, 肯定和N中第三行某列的元素相乘)
                int t;//和tp同样的套路
                if(brow < N.rows) t = N.rpos[brow+1];
                else t = N.elems+1;

                for(int q = N.rpos[brow]; q < t; q++){//这里不是直接算出M的某行乘N的某列的值
                    //而是用M的某行,去遍历N的所有行,然后对每行的应该在同列的值进行累加
                    int ccol = N.data[q].j;
                    ctemp[ccol] += M.data[p].e * N.data[q].e;//累加器
                }
            }//求得Q中第crow(=arow)行的非零元

            for(int ccol = 1; ccol <= Q.cols; ++ccol){//用M的一行遍历完了整个N矩阵
                if(ctemp[ccol]){
                    Q.data[++Q.elems] = {arow, ccol, ctemp[ccol]};//压缩存储
                }
            }
        }
    }
    return Q;
}

RLSMatrix makeMat(){
    /*
     * 输入rows, cols, 三元组个数
     * 输入三元组
     */
    RLSMatrix A;
    cin>>A.rows>>A.cols>>A.elems;
    for(int i = 1;i <= A.elems;i++){
        int x,y,e;
        cin>>x>>y>>e;
        A.data[i] = {x,y,e};
    }

    /*
     根据乘法中这段代码写的初始化rpos数组:
         int tp;      //用tp指向data中该行行末的下一个位置,方便遍历
         if(arow < M.rows) tp = M.rpos[arow+1];    //如果不是最后一行, tp指向下一行行首
         else tp = M.elems+1;   //最后一行, tp直接指向data中最后一个+1
         for(int p = M.rpos[arow];p < tp; p++) ...
        可以看出如果某行没有元素,则让tp = M.rpos[arow]则可不进行遍历,也就是M.rpos[arow] = M.rpos[arow+1]
        而当最后几行为空时,并不会执行else tp = M.elems+1;这时应该把rpos最后几行空的填成M.elems+1
     */

    int row = 1;
    for(int e = 1;e <= A.elems; e++){
        int arow = A.data[e].i;
        while(row <= arow){
            A.rpos[row++] = e;
        }
    }
    while(row <= A.rows){//如果最后几行没有元素,让索引指到elems+1的位置
        A.rpos[row++] = A.elems+1;
    }
    return A;
}

void printMat(RLSMatrix A){
    cout<<"rows:"<<A.rows<<" cols:"<<A.cols<<" elems:"<<A.elems<<endl;
    cout<<"-------------------------"<<endl;
    cout<<"data:"<<endl;
    for(int i = 1;i <= A.elems;i++) {
        cout<<setw(4)<<A.data[i].i<<setw(4)<<A.data[i].j<<setw(4)<<A.data[i].e<<endl;
    }
    cout<<"-------------------------"<<endl;
    cout<<"rpos: ";
    for(int j = 1;j <= A.rows; j++){
        cout<<A.rpos[j]<<' ';
    }
    cout<<endl<<endl;
}

int main(){
    ios::sync_with_stdio(false);
    RLSMatrix A = makeMat();
    printMat(A);
    RLSMatrix B = makeMat();
    printMat(B);
    RLSMatrix C = mul(A, B);
    printMat(C);
    return 0;
}
/*
 严书上的例子
 3 4 4
 1 1 3
 1 4 5
 2 2 -1
 3 1 2

 4 2 4
 1 2 2
 2 1 1
 3 1 -2
 3 2 4
 */

 

 
假设我们有两个矩阵A和B，我们必须找到AB的结果。我们可以假设A的列号等于B的行号。
 
因此，如果输入像[[1,0,0]，[-1,0,3]] [[7,0,0]，[0,0,0]，[0,0,1]] ，100
 
-103
 
700
 
000
 
001
 
那么输出将是[[7,0,0]，[-7,0,3]]700
 
-703
 
为了解决这个问题，我们将遵循以下步骤-r1：= A的大小，r2：= B的大小
 
c1：= A [0]的大小，c2：= B [0]的大小
 
定义顺序为r1 x c2的一个2D数组ret
 
定义成对的数组sparseA [r1]
 
对于初始化i：= 0，当i 
 
在sparseA [i]的末尾插入{j，A [i，j]}
 
对于初始化j：= 0，当j 
 
对于初始化i：= 0，当i 
 
ret [i，k]：= ret [i，k] + sparseA [i，j] * B [x，k]的第二个元素
 
x：= sparseA [i，j]的第一个元素
 
如果B [x，k]不等于0，则-
 
对于初始化j：= 0，当j 
 
返回ret
 
例
 
让我们看下面的实现以更好地理解-class Solution {
 
public:
 
vector multiply(vector& A, vector& B) {
 
int r1 = A.size();
 
int r2 = B.size();
 
int c1 = A[0].size();
 
int c2 = B[0].size();
 
vector  ret(r1, vector 
 
vector  > sparseA[r1];
 
for(int i = 0; i 
 
for(int j = 0; j 
 
if(A[i][j] != 0)sparseA[i].push_back({j, A[i][j]});
 
}
 
}
 
for(int i = 0; i 
 
for(int j = 0; j 
 
for(int k = 0; k 
 
int x = sparseA[i][j].first;
 
if(B[x][k] != 0){
 
ret[i][k] += sparseA[i][j].second * B[x][k];
 
}
 
}
 
}
 
}
 
return ret;
 
}
 
};
 
输入值{{1,0,0},{-1,0,3}},{{7,0,0},{0,0,0},{0,0,1}}
 
输出结果[[7, 0, 0, ],[-7, 0, 3, ],]

【数据结构】稀疏矩阵乘法
 
1.传统矩阵相乘的算法使用三个嵌套循环实现，算法复杂度为O（m * n1 * n2）
 2.使用三元组顺序表存储稀疏矩阵时，实现 Q= M * N，对于M中M（i，j）元素来说，只需要与N中第j行元素N（j,q）相乘，再存入Q（i，q）中。为了实现这一操作，增加一个向量rpos，表示每一行的第一个非零元在三元组中的位置，rpos作用相当于快速转置中的cpot向量。
 这种结构叫做 行链接的顺序表
 
typedef struct {
	Triple data[MAXSIZE + 1]; //非零三元组表，data0未用
	int rpos[MAXRC + 1];	  //各行第一个非零元的位置表
	int mu, nu, tu;			  //三元组行数，列数，非零元素值
}RLSMatrix;
 
 
在创建矩阵时，需要求得矩阵的rpos向量，求法与cpot一致。
 
Status CreateRLSMatrix(RLSMatrix &M) {
	cout << "输入矩阵行数，列数，非零元素个数：" << endl;
	cin >> M.mu >> M.nu >> M.tu;
	cout << "依次输入" << M.tu << "个元素的行数，列数，元素值:" << endl;
	for (int i = 1; i <= M.tu; i++)
	{
		cin >> M.data[i].i >> M.data[i].j >> M.data[i].e;
	}
	int* num = (int*)malloc(M.mu * sizeof(int));
	int arow;
	int q;
	if (M.tu) {
		for ( arow = 1; arow <= M.mu; ++arow) num[arow] = 0;
		for (int t = 1; t <= M.tu; ++t) ++num[M.data[t].i];//求M中每一行含非零元个数
		M.rpos[1] = 1;
		//求第col列中第一个非零元在b.data 中的序号
		for (arow = 2; arow <= M.mu; ++arow)M.rpos[arow] = M.rpos[arow - 1] + num[arow - 1];
	}
		cout << "矩阵构造完成" << endl;
	return OK;
}
 
 
稀疏矩阵乘法
 算法思路：
 大循环是对M中的每一行进行。
 rpos向量用来确定每行元素个数
 再对该行每个元素M（i，j）进行处理，找到N中对应的第 j 行，将这一行所有元素与M（i，j）相乘，结果存入ctemp累加器中，该行元素都处理完成后，ctemp中所存储的便是Q中第 i 行数据
 
Status MultSMatrix(RLSMatrix M, RLSMatrix N, RLSMatrix &Q) {
	//求稀疏矩阵乘积 Q = M x N  采用行逻辑链接存储表示
	if (M.nu != N.mu)return ERROR;
	Q.mu = M.mu; Q.nu = N.nu; Q.tu = 0;
	int arow,tp,p,brow,t,q,ccol;
	int* ctemp = (int*)malloc((M.mu+1) * sizeof(int));
	if (M.tu * N.tu != 0) {
		for (arow = 1; arow <= M.mu; ++arow) {
		//对M的每行进行处理
			for (int i = 0; i <= M.mu; i++)ctemp[i] = 0;//累加器清零
			Q.rpos[arow] = Q.tu + 1;
			if (arow < M.mu)tp = M.rpos[arow + 1];
			else tp = M.tu + 1;//确定M中该行元素个数
			for ( p = M.rpos[arow]; p < tp; p++)
			{//M中对该行的每个元素进行处理
				brow = M.data[p].j;
				if (brow < N.mu) t = N.rpos[brow + 1];
				else t = N.tu + 1;
				for (q = N.rpos[brow]; q < t; ++q) {
				//将该元素与矩阵N中对应行的元素进行相乘，结果存入累加器
					ccol = N.data[q].j;
					ctemp[ccol] += M.data[p].e *N.data[q].e;
				}
			}
			for(ccol =1;ccol<=Q.nu;++ccol)
				if (ctemp[ccol]) {
				//第arow行处理完毕，将累加器中值存入Q中
					if (++Q.tu > MAXSIZE)return ERROR;
					Q.data[Q.tu] = { arow, ccol, ctemp[ccol] };
				}//if
		}//for arow
	}//if
}
 
 
该算法总的时间复杂度为O（M.mu * N.nu + M.tu * N.tu/N.mu）
 ctemp初始化 O（M.mu * N.nu）
 求Q中所有非零元 O（M.tu * N.tu / N.mu）
 压缩存储 O（M.mu * N.nu）