线性矩阵方程的最小二乘解返回函数

Question

我是一名有用的助手，可以为您翻译文本。

我一直在尝试将代码从Python重写为Swift，但我卡在了一个函数上，该函数应返回线性矩阵方程的最小二乘解。 有人知道是否有一个用Swift编写的库，其中具有与{{link1：numpy.linalg.lstsq}}等效的方法吗？ 我将非常感激您的帮助。

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代码：

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]

- wtznc

---

添加代码在这里！没有代码无法帮助。 - Mtoklitz113

加速框架包括 BLAS 库，其中有用于线性最小二乘问题的函数。然而，要从 Swift 中使用这些函数需要一些工作 :) - Martin R

很遗憾，目前还没有一种方法可以解决LLS问题。 - wtznc

难道不是 DGELS（或其变体）是您所需的吗？ - Martin R

你能提及一下用于计算线性矩阵方程的最小二乘解的公式吗？ - Mtoklitz113

有一些注意事项，但如果您在2个变量中具有超定的n个方程组，即如果您的问题可以写成A * x = b，其中A的形状为(n, 2)，x的形状为(2,)，b的形状为(n,)，则最小二乘解与等效的A.T * A * x = A.T * b的解相同，由于A.T * A的形状为(2, 2)，而A.T * b的形状为(2,)，因此您有一个非常简单的系统，可以轻松解决，例如使用Cramer's rule。从头开始编写应该非常简单。 - Jaime

Answer 1

Accelerate框架包含LAPACK线性代数包，其中包括DGELS函数，用于解决欠定或超定的线性系统。根据文档：

DGELS通过使用A的QR或LQ分解来解决涉及M×N矩阵A或其转置的过定或欠定实线性系统。假设A具有完全秩。

以下是如何从Swift中使用该函数的示例。这本质上是此C示例代码的翻译。

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：欠定系统，对于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)))
}

- Martin R

---

这太棒了！谢谢你！ - gtmtg