¿Cómo resolver un sistema lineal de matrices en Scala Breeze?

Question

¿Cómo resolver un sistema lineal de matrices en Scala Breeze? es decir, tengo Ax = b, donde A es una matriz (generalmente positiva definida), y x y b son vectores.

Veo que hay una descomposición de Cholesky disponible, pero parece que no puedo encontrar un solucionador? (si fuera matlab podría hacer x = b \ A. Si fuera poco fácil podría hacer x = A.solve (b))

Answer 1

Al parecer, es bastante simple en realidad, y construido en Scala-brisa como un operador:

x = A \ b 

Eso no utilizar Cholesky, que utiliza la descomposición LU, que es creo que alrededor de la mitad tan rápido, pero ambos son O (n^3), tan comparables.

Answer 2

Bueno, al final escribí mi propio solucionador. No estoy seguro de si esta es la forma óptima de hacerlo, pero no parece irrazonable. :

// 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 
    } 
}