稀疏矩阵与稠密矩阵以及 LibRec 中的实现

稀疏矩阵

在矩阵中，若数值为0的元素数目远远多于非0元素的数目时，则称该矩阵为稀疏矩阵。与之相反，若非0元素数目占大多数时，则称该矩阵为稠密矩阵。

注：压缩存储的矩阵可以分为特殊矩阵和稀疏矩阵。对于那些具有相同元素或零元素在矩阵中分布具有一定规律的矩阵，被称之为特殊矩阵。对于那些零元素数据远远多于非零元素数目，并且非零元素的分布没有规律的矩阵称之为稀疏矩阵。

稀疏矩阵的特性

稀疏矩阵其非零元素的个数远远小于零元素的个数，而且这些非零元素的分布也没有规律。
稀疏因子是用于描述稀疏矩阵的非零元素的比例情况。设一个 n*m 的稀疏矩阵 A 中有 t 个非零元素，则稀疏因子δ 的计算公式：$δ=t/mn$ (当这个值小于等于0.05时，可以认为是稀疏矩阵)*

稀疏矩阵的压缩存储

存储矩阵的一般方法是采用二维数组，其优点是可以随机地访问每一个元素，因而能够较容易地实现矩阵的各种运算，如转置运算、加法运算、乘法运算等。但对于稀疏矩阵来说，采用二维数组的存储方法既浪费大量的存储单元用来存放零元素，又要在运算中花费大量的时间来进行零元素的无效计算，显然不科学。所以必须考虑对稀疏矩阵进行压缩存储。

对稀疏矩阵进行压缩存储的一种较好的方法是：只存储在矩阵中极少数的非零元素，为此必须对每一个非零元素，保存它的下标和值。可以采用一个三元组Trituple<row,column,value>来唯一地确定一个矩阵元素。因此，稀疏矩阵需要使用一个三元组数组(亦称为三元组表)来表示。

常见的压缩存储方式有压缩行和压缩列。

对于如下矩阵

以三元组的方式(假设索引从1开始)为：

(row_index,col_index,value)

按行存储
[(1,1,10),(1,5,-2),			//第1行
(2,1,3), (2,2,9),(2,6,3),	//第2行
(3,2,7),(3,3,8),(3,4,7),	//第3行
(4,1,3),(4,5,5),			//第4行
(5,2,8),(5,4,9),(5,5,9),	//第5行
(6,2,4),(6,5,2),(6,6,-1)]	//第6行
 	 
按列存储
[(1,1,10),(2,1,3),(4,1,3),			//第1列
(2,2,9),(3,2,7),(5,2,8),(6,2,4),	//第2列
(3,3,8),							//第3列
(3,4,7),(5,4,9),					//第4列
(1,5,-2),(4,5,5),(5,5,9),(6,5,2),	//第5列
(2,6,3),(6,6,-1)]					//第6列

Compressed Row Storage (CRS) 压缩行存储

压缩行和列存储格式是最通用的：它们不对矩阵的稀疏性结构做出假设，并且不存储任何不必要的元素。但是，效率不高，需要对矩阵向量积或预处理器求解中的每个标量运算进行间接寻址步骤。

CSR矩阵格式通过用三个一维的数组来存储一个m×n的矩阵M。定义nnz（Num-non-zero）为矩阵M中非0元素的个数。

val[] : 存储非 0 元素，长度为 nnz。
col_ind[]：存储 val 中元素的列索引，nnz。
row_ptr[]：存储每一行中第一个非零元素在 val 中的下标，长度为 m+1,其中 $row_ptr[m+1] = nnz+1$ 。这里即将行元素索引压缩了。

我们只需要 $2*nnz+m+1$的空间去存储矩阵。

index     = [1,  2, 3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14, 15,16] //数组下标
value     = [10,-2, 3, 9, 3, 7, 8, 7, 3, 5, 8, 9, 9, 4, 2, -1] //按行读取时，依次出现的非零值
col_index = [1,  5, 1, 2, 6, 2, 3, 4, 1, 5,  2, 4, 5, 2, 6,  6] //按行读取时，依次出现的非零值所在的列下标
row_ptr   = [1,  3, 6, 9,11,14,17] //按行读取时，每一行中第一个出现的非零值在value数组中的下标值,可以参考三元组中按行存储。

Compressed Column Storage (CCS) 压缩列存储

index     = [1,  2, 3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14, 15,16] //数组下标
value     = [10, 3, 3, 9, 7, 8, 4, 8, 7, 9, -2, 5, 9, 2, 3, -1] //按列读取时，依次出现的非零值
row_index = [1,  2, 4, 2, 3, 5, 6, 3, 3, 5,  1, 4, 5, 6, 2,  6] //案列读取时，依次出现非零值所在的行下标
col_ptr   = [1, 4, 8, 9,11,15,17] //按列读取时，每一列中第一个出现的非零值在value数组中的下标值,可以参考三元组中按列存储。

LibRec 实现

由于推荐系统场景中，用户物品的行为矩阵是一个稀疏矩阵，LibRec 2.0 中的 SparseMatrix.java类实现了这一数据结构。其同时实现了 CRS 与 CCS 压缩存储。

我们先查看 LibRec 如何读取输入的用户评分/行为矩阵。以 TextDataConvertor.java为例，readData()函数实现了数据读取，存入 SparseMatrx。主要查看以下部分。

// Table {row-id, col-id, rate}  矩阵的三元组
Table<Integer, Integer, Double> dataTable = HashBasedTable.create();
// Map {col-id, multiple row-id}: used to fast build a rating matrix 是一个Multimap，存储列索引对应的行索引用以加快稀疏矩阵构建
Multimap<Integer, Integer> colMap = HashMultimap.create();
....
    // 不断读取数据行
		for (int i = 0; i < loopLength; i++) {
                    String line = bufferData[i];
                    String[] data = line.trim().split("[ \t,]+");
                    String user = data[0];
                    String item = data[1];
                    Double rate = ((dataColumnFormat.equals("UIR") || 		dataColumnFormat.equals("UIRT")) && data.length >= 3) ? Double.valueOf(data[2]) : 1.0;
                    // binarize the rating for item recommendation task
                    if (binThold >= 0) {
                        rate = rate > binThold ? 1.0 : 0.0;
                    }
                    // inner id starting from 0 
            		//稀疏矩阵的行索引，以判断 userId 是否存在来递增
                    int row = userIds.containsKey(user) ? userIds.get(user) : userIds.size();
                    userIds.put(user, row);
				  //稀疏矩阵的列索引，以判断 itemId 是否存在来递增
                    int col = itemIds.containsKey(item) ? itemIds.get(item) : itemIds.size();
                    itemIds.put(item, col);
                    dataTable.put(row, col, rate);
                    colMap.put(col, row);
                    // record rating's issuing time
                    if (StringUtils.equals(dataColumnFormat, "UIRT") && data.length >= 4) {
                        if (timeTable == null) {
                            timeTable = HashBasedTable.create();
                        }
                        // convert to million-seconds
                        long mms = 0L;
                        try {
                            mms = Long.parseLong(data[3]); // cannot format
                            // 9.7323480e+008
                        } catch (NumberFormatException e) {
                            mms = (long) Double.parseDouble(data[3]);
                        }
                        long timestamp = timeUnit.toMillis(mms);
                        timeTable.put(row, col, timestamp);
                    }
                }
                if (!isComplete) {
                    bufferLine = bufferData[bufferData.length - 1];
                } else {
                    bufferLine = "";
                }
                buffer.clear();
            }
            fileRead.close();
            fis.close();
        }
        int numRows = numUsers(), numCols = numItems();
        // build rating matrix  
        preferenceMatrix = new SparseMatrix(numRows, numCols, dataTable, colMap);
        if (timeTable != null)
            datetimeMatrix = new SparseMatrix(numRows, numCols, timeTable, colMap);
        // release memory of data table
        dataTable = null;
        timeTable = null;
    }

我们通过构造器来查看如何进行稀疏矩阵构建。

/**
     * Construct a sparse matrix with both CRS and CCS structures
     *
     * @param rows       number of rows
     * @param cols       number of columns
     * @param dataTable  data table  传入原始数据矩阵的三元组
     * @param colMap     column map  存储列索引对应的行索引用以加快稀疏矩阵构建
     */
    public SparseMatrix(int rows, int cols, Table<Integer, Integer, ? extends Number> dataTable, Multimap<Integer, Integer> colMap) {
        numRows = rows;
        numColumns = cols;
        valueSet = new TreeSet<>();
        construct(dataTable, colMap);
    }
 /**
     * Construct a sparse matrix
     *
     * @param dataTable       data table
     * @param columnStructure column structure
     */
    private void construct(Table<Integer, Integer, ? extends Number> dataTable,
                           Multimap<Integer, Integer> columnStructure) {
        int nnz = dataTable.size();
        // CRS 压缩行存储
        rowPtr = new int[numRows + 1];
        colInd = new int[nnz];
        rowData = new double[nnz];
        int j = 0;
        for (int i = 1; i <= numRows; ++i) {
            // 存储 rowPtr，为稀疏矩阵中每行第一个非0元素在 rowData 中的索引
            Set<Integer> cols = dataTable.row(i - 1).keySet();
            rowPtr[i] = rowPtr[i - 1] + cols.size();
            // 存储 colInd，按行读取时，依次出现的非零值所在的列下标
            for (int col : cols) {
                colInd[j++] = col;
                if (col < 0 || col >= numColumns)
                    throw new IllegalArgumentException("colInd[" + j + "]=" + col
                            + ", which is not a valid column index");
            }
            // 排序，保证列是从小到大的顺序
            Arrays.sort(colInd, rowPtr[i - 1], rowPtr[i]);
        }
        // CCS
        colPtr = new int[numColumns + 1];
        rowInd = new int[nnz];
        colData = new double[nnz];
        j = 0;
        for (int i = 1; i <= numColumns; ++i) {
            // dataTable.col(i-1) is more time-consuming than columnStructure.get(i-1)
            // 由于 dataTable.col(i-1) 比 columnStructure.get(i-1) 慢很多，使用 colMap 来加快构建速度
            Collection<Integer> rows = columnStructure != null ? columnStructure.get(i - 1) : dataTable.column(i - 1)
                    .keySet();
            colPtr[i] = colPtr[i - 1] + rows.size();
            for (int row : rows) {
                rowInd[j++] = row;
                if (row < 0 || row >= numRows)
                    throw new IllegalArgumentException("rowInd[" + j + "]=" + row + ", which is not a valid row index");
            }
            Arrays.sort(rowInd, colPtr[i - 1], colPtr[i]);
        }
        // set data  设置稀疏矩阵中的非0值，包括 rowData 和 colData
        for (Cell<Integer, Integer, ? extends Number> en : dataTable.cellSet()) {
            int row = en.getRowKey();
            int col = en.getColumnKey();
            double val = en.getValue().doubleValue();
            set(row, col, val);
        }
    }

以设置矩阵非0值与通过行列索引查找 rowData 的相应索引位置

/**
     * Set a value to entry [row, column]
     *
     * @param row    row id
     * @param column column id
     * @param val    value to set
     */
    public void set(int row, int column, double val) {
        int index = getCRSIndex(row, column);
        rowData[index] = val;
        index = getCCSIndex(row, column);
        colData[index] = val;
        valueSetAdd(val);
    }
/**
     * Finds the insertion index of CRS
     * 通过行列索引查找 rowData 的相应索引位置
     * @param row  the index of row
     * @param col  the index of column
     */
    private int getCRSIndex(int row, int col) {
        // 二分查找， colInd 存储 rowData 对应的矩阵列下标，rowPtr[row] 代表了 colIndex 中第 row 行第一个非零元素对应的位置，col 是需要查找的值
        int i = Arrays.binarySearch(colInd, rowPtr[row], rowPtr[row + 1], col);
        if (i >= 0 && colInd[i] == col)
            return i;
        else
            throw new IndexOutOfBoundsException("Entry (" + (row + 1) + ", " + (col + 1)
                    + ") is not in the matrix structure");
    }

查找稀疏矩阵中的对应元素

/**
     * Retrieve value at entry [row, column]
     *
     * @param row    row id
     * @param column column id
     * @return value at entry [row, column]
     */
    public double get(int row, int column) {
		// 与 getCRSIndex 原理相同，通过 rowPtr 和 colInd 来进行 rowData 的元素定位
        int index = Arrays.binarySearch(colInd, rowPtr[row], rowPtr[row + 1], column);
        if (index >= 0)
            return rowData[index];
        else
            return 0;
    }

参考

Compressed Row Storage (CRS)

Compressed Column Storage (CCS)

LibRec

http://westerly-lzh.github.io/cn/2014/12/Sparse-Matrix-Storage-Style/