2012-09-28 14 views
6

Scala esintisinde lineer bir matris sistemi nasıl çözülür? yani, A = matrisi (genellikle pozitif bir kesin) olan ve x ve b vektörleri olan Ax = b var.Scala esintisinde lineer bir matris sistemi nasıl çözülür?

Mevcut bir ayrışma olduğunu görebiliyorum ama çözücü bulamıyorum. (eğer matlab olsaydı, x = b \ A yapabilirdim. Eğer bu bir scipy ise x = A.solve (b))

cevap

6

Görünüşe göre, aslında oldukça basittir ve bir operatör olarak Scala-esinti yerleşik:

x = A \ b 

O Cholesky kullanmak etmez, o LU bozulmasını kullanır, ben yarım düşünmek olduğunu hızlı, ama ikisi de O (n^3), yani karşılaştırılabilir.

4

Sonunda kendi çözücüyü yazdım. Bunun yapmanın en iyi yolu olup olmadığından emin değilim, ama mantıksız görünmüyor mu? :

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