scalalinear-algebrascalala

How to do X * diag(Y) in Scala Breeze?


How to do X * diag(Y) in Scala Breeze? X could be for example a CSCMatrix and Y could be a DenseVector?

In MATLAB syntax, this would be:

X * spdiags(0, Y, N, N )

Or:

X .* repmat( Y', K, 0 )

In SciPy syntax, this would be a 'broadcast multiply':

Y * X

How to do X * diag(Y) in Scala Breeze?


Solution

  • I wrote my own sparse diagonal method, and dense / sparse multiplication method in the end.

    Use like this:

    val N = 100000
    val K = 100
    val A = DenseMatrix.rand(N,K)
    val b = DenseVector.rand(N)
    val c = MatrixHelper.spdiag(b)
    val d = MatrixHelper.mul( A.t, c )
    

    Here are the implementations of spdiag and mul:

    // Copyright Hugh Perkins 2012
    // You can use this under the terms of the Apache Public License 2.0
    // http://www.apache.org/licenses/LICENSE-2.0
    
    package root
    
    import breeze.linalg._
    
    object MatrixHelper {
       // it's only efficient to put the sparse matrix on the right hand side, since 
       // it is a column-sparse matrix
       def mul( A: DenseMatrix[Double], B: CSCMatrix[Double] ) : DenseMatrix[Double] = {
          val resultRows = A.rows
          val resultCols = B.cols
          var row = 0
          val result = DenseMatrix.zeros[Double](resultRows, resultCols )
          while( row < resultRows ) {
             var col = 0
             while( col < resultCols ) {
                val rightRowStartIndex = B.colPtrs(col)
                val rightRowEndIndex = B.colPtrs(col + 1) - 1
                val numRightRows = rightRowEndIndex - rightRowStartIndex + 1
                var ri = 0
                var sum = 0.
                while( ri < numRightRows ) {
                   val inner = B.rowIndices(rightRowStartIndex + ri)
                   val rightValue = B.data(rightRowStartIndex + ri)
                   sum += A(row,inner) * rightValue
                   ri += 1
                }
                result(row,col) = sum
                col += 1
             }
             row += 1
          }
          result
       }   
    
       def spdiag( a: Tensor[Int,Double] ) : CSCMatrix[Double] = {
          val size = a.size
          val result = CSCMatrix.zeros[Double](size,size)
          result.reserve(a.size)
          var i = 0
          while( i < size ) {
             result.rowIndices(i) = i
             result.colPtrs(i) = i
             result.data(i) = i
             //result(i,i) = a(i)
             i += 1
          }
          //result.activeSize = size
          result.colPtrs(i) = i
          result
       }
    }