加速框架包括LAPACK http://www.netlib.org/lapack/index.html线性代数包,
其中有一个DGELS http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga1df516c81d3e902cca1fc79a7220b9cb.html#ga1df516c81d3e902cca1fc79a7220b9cb函数来求解欠定或超定线性系统。从文档中:
DGELS 求解超定或欠定实线性系统
涉及 M×N 矩阵 A 或其转置,使用 QR 或 LQ
A 的因式分解。假设 A 具有满秩。
以下是如何从 Swift 使用该函数的示例。
它本质上是一个翻译这个C示例代码 https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgels_ex.c.htm.
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
and CInt
必要的。
示例1:超定系统,由http://www.seas.ucla.edu/~vandenbe/103/lectures/ls.pdf 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 and 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)))
}