Làm thế nào để giải quyết một hệ thống tuyến tính của ma trận trong gió scala?

Question

Làm thế nào để giải quyết một hệ thống tuyến tính của ma trận trong gió scala? tức là, tôi có Ax = b, trong đó A là ma trận (thường là xác định dương), và x và b là vectơ.

Tôi có thể thấy rằng có một phân tích có sẵn, nhưng tôi dường như không thể tìm thấy người giải quyết? (Nếu nó là matlab tôi có thể làm x = b \ A. Nếu nó đã được scipy tôi có thể làm x = A.solve (b))

Answer 1

Rõ ràng, nó là khá đơn giản trong thực tế, và được xây dựng vào scala-gió như một nhà điều hành:

x = A \ b 

Nó không sử dụng Cholesky, nó sử dụng phân tích lu, đó là tôi nghĩ về một nửa nhanh, nhưng chúng đều là O (n^3), nên có thể so sánh được.

Answer 2

Vâng, tôi đã viết giải quyết của riêng tôi cuối cùng. Tôi không chắc đây có phải là cách tối ưu để làm điều đó, nhưng nó có vẻ không hợp lý? :

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