wtznc - 2 years ago 128
Swift Question

# Function which returns the least-squares solution to a linear matrix equation

I must rewrite code from Python to Swift but I'm stuck on function which should return the least-squares solution to a linear matrix equation. Does anyone of you know a library written in Swift which has an equivalent method to the

? I'd be grateful for your help.

Python code:

``````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]
``````

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 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.

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)))
}
``````

Some notes:

• `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 `var`s.
• A C integer is a 32-bit integer, which makes some conversions between `Int` and `CInt` necessary.

Example 1: Overdetermined system, from 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]
}
``````

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]
}
``````
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download