How to solve a linear system of matrices in scala breeze?

Question

How to solve a linear system of matrices in scala breeze? ie, I have Ax = b, where A is a matrix (usually positive definite), and x and b are vectors.

I can see that there is a cholesky decomposition available, but I couldn't seem to find a solver? (if it was matlab I could do x = b \ A. If it was scipy I could do x = A.solve(b) )

Legitimist

Answer 1

Apparently, it is quite simple in fact, and built into scala-breeze as an operator:

x = A \ b

It doesnt use Cholesky, it uses LU decomposition, which is I think about half as fast, but they are both O(n^3), so comparable.

Legitimist

Answer 2

Well, I wrote my own solver in the end. I'm not sure if this is the optimal way to do it, but it doesn't seem unreasonable? :

// 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 Solver {
   // solve Ax = b, for x, where A = choleskyMatrix * choleskyMatrix.t
   // choleskyMatrix should be lower triangular
   def solve( choleskyMatrix: DenseMatrix[Double], b: DenseVector[Double] ) : DenseVector[Double] = {
      val C = choleskyMatrix
      val size = C.rows
      if( C.rows != C.cols ) {
          // throw exception or something
      }
      if( b.length != size ) {
          // throw exception or something
      }
      // first we solve C * y = b
      // (then we will solve C.t * x = y)
      val y = DenseVector.zeros[Double](size)
      // now we just work our way down from the top of the lower triangular matrix
      for( i <- 0 until size ) {
         var sum = 0.
         for( j <- 0 until i ) {
            sum += C(i,j) * y(j)
         }
         y(i) = ( b(i) - sum ) / C(i,i)
      }
      // now calculate x
      val x = DenseVector.zeros[Double](size)
      val Ct = C.t
      // work up from bottom this time
      for( i <- size -1 to 0 by -1 ) {
         var sum = 0.
         for( j <- i + 1 until size ) {
            sum += Ct(i,j) * x(j)
         }
         x(i) = ( y(i) - sum ) / Ct(i,i)
      }
      x
   }
}

Legitimist

In case you really need to solve via a Cholesky decomposition. Breeze sits on top of LAPACK. There are two subroutines there for dealing with the Cholesky decomposition. ?potrf (with ? being s for single-precision, d for double-precision, or c for complex) does the decomposition. ?potrs uses the decomposition to actually solve the equation. Like cholesky() in Breeze is a thin wrapper over dpotrf, you can simply call dpotrs instead of writing your own solver. – Levi