최소 제곱 해를 선형 행렬 방정식으로 반환하는 함수

Question

파이썬에서 스위프트로 코드를 다시 작성해야하지만 선형 행렬 방정식에 대한 최소 제곱 해를 반환해야하는 함수가 필요합니다. Swift에 쓰여진 도서관을 아는 사람이 누구입니까? numpy.linalg.lstsq과 같은 방법이 있습니까? 당신의 도움에 감사드립니다.

파이썬 코드 : 지금까지

a = numpy.array([[p2.x-p1.x,p2.y-p1.y],[p4.x-p3.x,p4.y-p3.y],[p4.x-p2.x,p4.y-p2.y],[p3.x-p1.x,p3.y-p1.y]]) 
b = numpy.array([number1,number2,number3,number4]) 
res = numpy.linalg.lstsq(a,b) 
result = [float(res[0][0]),float(res[0][1])] 
return result 

스위프트 코드 :

var matrix1 = [[p2.x-p1.x, p2.y-p1.y],[p4.x-p3.x, p4.y-p3.y], [p4.x-p2.x, p4.y-p2.y], [p3.x-p1.x, p3.y-p1.y]] 
var matrix2 = [number1, number2, number3, number4] 

Answer 1

는 프레임 워크가 과소 또는 중복 - 결정된 선형 시스템을 해결하기 위해 DGELS 기능을 가지고 LAPACK 선형 대수 패키지, 을 포함 가속화 . 문서에서 :

DGELS 전체 것으로 가정하는 QR 또는 A.의 LQ 인수 분해를 이용하여 M-바이 N 행렬 A, 또는 그 전치를 포함하는 중복 - 결정된 또는 underdetermined 실제 선형 시스템 해결 계급.

다음은 Swift에서이 기능을 어떻게 사용할 수 있는지 예입니다. 이것은 본질적으로 this C sample code의 번역입니다.

func solveLeastSquare(A A: [[Double]], B: [Double]) -> [Double]? { 
    precondition(A.count == B.count, "Non-matching dimensions") 

    var mode = Int8(bitPattern: UInt8(ascii: "N")) // "Normal" mode 
    var nrows = CInt(A.count) 
    var ncols = CInt(A[0].count) 
    var nrhs = CInt(1) 
    var ldb = max(nrows, ncols) 

    // Flattened columns of matrix A 
    var localA = (0 ..< nrows * ncols).map { 
     A[Int($0 % nrows)][Int($0/nrows)] 
    } 

    // Vector B, expanded by zeros if ncols > nrows 
    var localB = B 
    if ldb > nrows { 
     localB.appendContentsOf([Double](count: ldb - nrows, repeatedValue: 0.0)) 
    } 

    var wkopt = 0.0 
    var lwork: CInt = -1 
    var info: CInt = 0 

    // First call to determine optimal workspace size 
    dgels_(&mode, &nrows, &ncols, &nrhs, &localA, &nrows, &localB, &ldb, &wkopt, &lwork, &info) 
    lwork = Int32(wkopt) 

    // Allocate workspace and do actual calculation 
    var work = [Double](count: Int(lwork), repeatedValue: 0.0) 
    dgels_(&mode, &nrows, &ncols, &nrhs, &localA, &nrows, &localB, &ldb, &work, &lwork, &info) 

    if info != 0 { 
     print("A does not have full rank; the least squares solution could not be computed.") 
     return nil 
    } 
    return Array(localB.prefix(Int(ncols))) 
} 

일부 메모 :

• dgels_()는 전달 매트릭스 및 벡터 데이터를 수정하고 에게 A의 열을 포함하는 "평면"배열로 행렬을 기대하고있다. 또한 오른쪽은 길이가 max(M, N) 인 배열로 예상됩니다. 이러한 이유로 입력 데이터는 먼저 로컬 변수에 복사됩니다.
• 모든 인수는 dgels_()을 참조해야하며, 따라서 은 모두 var에 저장됩니다.
• C 정수는 32 비트 정수로 Int과 CInt 사이의 변환이 필요합니다.

예 1에서 http://www.seas.ucla.edu/~vandenbe/103/lectures/ls.pdf 중복 - 결정된 시스템.

let A = [[ 2.0, 0.0 ], 
     [ -1.0, 1.0 ], 
     [ 0.0, 2.0 ]] 
let B = [ 1.0, 0.0, -1.0 ] 
if let x = solveLeastSquare(A: A, B: B) { 
    print(x) // [0.33333333333333326, -0.33333333333333343] 
} 

예 2 Underdetermined 시스템 최소 규범에 x_1 + x_2 + x_3 = 1.0 용액.

let A = [[ 1.0, 1.0, 1.0 ]] 
let B = [ 1.0 ] 
if let x = solveLeastSquare(A: A, B: B) { 
    print(x) // [0.33333333333333337, 0.33333333333333337, 0.33333333333333337] 
} 

---

업데이트 스위프트 3 및 스위프트 4 :

func solveLeastSquare(A: [[Double]], B: [Double]) -> [Double]? { 
    precondition(A.count == B.count, "Non-matching dimensions") 

    var mode = Int8(bitPattern: UInt8(ascii: "N")) // "Normal" mode 
    var nrows = CInt(A.count) 
    var ncols = CInt(A[0].count) 
    var nrhs = CInt(1) 
    var ldb = max(nrows, ncols) 

    // Flattened columns of matrix A 
    var localA = (0 ..< nrows * ncols).map { (i) -> Double in 
     A[Int(i % nrows)][Int(i/nrows)] 
    } 

    // Vector B, expanded by zeros if ncols > nrows 
    var localB = B 
    if ldb > nrows { 
     localB.append(contentsOf: [Double](repeating: 0.0, count: Int(ldb - nrows))) 
    } 

    var wkopt = 0.0 
    var lwork: CInt = -1 
    var info: CInt = 0 

    // First call to determine optimal workspace size 
    var nrows_copy = nrows // Workaround for SE-0176 
    dgels_(&mode, &nrows, &ncols, &nrhs, &localA, &nrows_copy, &localB, &ldb, &wkopt, &lwork, &info) 
    lwork = Int32(wkopt) 

    // Allocate workspace and do actual calculation 
    var work = [Double](repeating: 0.0, count: Int(lwork)) 
    dgels_(&mode, &nrows, &ncols, &nrhs, &localA, &nrows_copy, &localB, &ldb, &work, &lwork, &info) 

    if info != 0 { 
     print("A does not have full rank; the least squares solution could not be computed.") 
     return nil 
    } 
    return Array(localB.prefix(Int(ncols))) 
}