2012-09-28 90 views
6

स्कैला हवा में मैट्रिस की रैखिक प्रणाली को कैसे हल करें? यानी, मेरे पास एक्स = बी है, जहां ए एक मैट्रिक्स (आमतौर पर सकारात्मक निश्चित) है, और एक्स और बी वैक्टर हैं।स्कैला हवा में मैट्रिस की रैखिक प्रणाली को कैसे हल करें?

मैं देख सकता हूं कि एक चापलूसी अपघटन उपलब्ध है, लेकिन मुझे एक सॉल्वर नहीं लग रहा था? (अगर यह मैटलैब था तो मैं x = b \ ए कर सकता था। अगर यह छोटा था तो मैं x = A.solve (बी) कर सकता था)

उत्तर

6

जाहिर है, यह वास्तव में बहुत सरल है, और एक ऑपरेटर के रूप में स्केला-हवा में बनाया गया:

x = A \ b 

यह Cholesky का उपयोग नहीं करता है, यह LU अपघटन का उपयोग करता है, जो कि मैं के बारे में आधे लगता है तेज़ी से, लेकिन वे ओ (एन^3) दोनों हैं, तो तुलनात्मक।

4

ठीक है, मैंने अंत में अपना स्वयं का सॉल्वर लिखा था। मुझे यकीन नहीं है कि यह करने का यह सबसे अच्छा तरीका है, लेकिन यह अनुचित नहीं लगता है? :

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