2012-09-29 32 views
0

如何在斯卡拉微風X * diag(Y)X可能是例如CSCMatrixY可能是DenseVector如何在Scala微風中做X * diag(Y)?

在MATLAB中的語法,這將是:

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

或者:

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

在SciPy的語法,這將是一個 '廣播乘法':

Y * X 

如何在Scala Breeze中做X * diag(Y)

回答

2

我寫了自己的稀疏對角線方法和最後的密/稀疏乘法方法。

使用這樣的:

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) 

這裏是SPDiag中和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 
    } 
}