将最小二乘解返回线性矩阵方程的函数 [英] Function which returns the least-squares solution to a linear matrix equation

问题描述

我一直在尝试将代码从Python重写为Swift，但我坚持使用该函数，该函数应将最小二乘解返回线性矩阵方程.有谁知道用Swift编写的库，它具有与

I have been trying to rewrite the code from Python to Swift but I'm stuck on the function which should return the least-squares solution to a linear matrix equation. Does anyone know a library written in Swift which has an equivalent method to the numpy.linalg.lstsq? I'd be grateful for your help.

Python代码:

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

到目前为止的快速代码:

Swift code so far:

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]

推荐答案

加速框架包括 LAPACK 线性代数包， 其中具有 Dellyer> -或超定线性系统.从文档中:

The Accelerate framework included the LAPACK linear algebra package, which has a DGELS function to solve under- or overdetermined linear systems. From the documentation:

DGELS解决过高或过低的实线性系统 使用QR或LQ涉及M×N矩阵A或其转置 A的因式分解.假定A具有最高等级.

DGELS solves overdetermined or underdetermined real linear systems involving an M-by-N matrix A, or its transpose, using a QR or LQ factorization of A. It is assumed that A has full rank.

下面是一个示例，说明如何从Swift中使用该功能. 它实质上是此C的翻译.示例代码.

Here is an example how that function can be used from Swift. It is essentially a translation of 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 s中.
• C整数是32位整数，它在 Int和CInt必需.

• dgels_() modifies the passed matrix and vector data, and expects the matrix as "flat" array containing the columns of A. Also the right-hand side is expected as an array with length max(M, N). For this reason, the input data is copied to local variables first.
• All arguments must be passed by reference to dgels_(), that's why they are all stored in vars.
• A C integer is a 32-bit integer, which makes some conversions between Int and CInt necessary.

示例1::从 示例2:系统欠佳，最低标准 x_1 + x_2 + x_3 = 1.0的解决方案.

Example 2: Underdetermined system, minimum norm solution to 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]
}

---

更新 Swift 3 和 Swift 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)))
}

这篇关于将最小二乘解返回线性矩阵方程的函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！