implement svd

Closes #33

README.md

@@ -453,6 +453,33 @@ print(P * L * U)
 
 See also `luInPlace()` for more advanced use cases that avoid creating full matrices.
 
+### Singular Value Decomposition (SVD)
+
+```swift
+let A = Matrix<Double>(NdArray.range(from: 1, to: 9).reshaped([2, 4]))
+let (U, s, Vt) = try A.svd()
+print(U)
+// [[-0.3761682344281408, -0.9265513797988838],
+// [-0.9265513797988838, 0.3761682344281408]]
+print(s)
+// [14.227407412633742, 1.2573298353791098]
+print(Vt)
+// [[ -0.3520616924890126, -0.44362578258952023, -0.5351898726900277, -0.6267539627905352],
+// [ 0.7589812676751458, 0.32124159914593237, -0.1164980693832819, -0.554237737912496],
+// [ -0.4000874340557387, 0.25463292200666415, 0.6909964581538871, -0.5455419461048127],
+// [ -0.3740722458438949, 0.7969705609558909, -0.471724384380099, 0.04882606926810252]]
+let Sd = Matrix(diag: s)
+let S = Matrix<Double>.zeros(A.shape)
+let mn = A.shape.min()!
+S[..<mn, ..<mn] = Sd
+print(S)
+// [[14.227407412633742, 0.0, 0.0, 0.0],
+// [ 0.0, 1.2573298353791098, 0.0, 0.0]]
+print(U * S * Vt)
+// [[1.0000000000000004, 2.0, 3.0000000000000004, 3.9999999999999996],
+// [ 4.999999999999999, 6.000000000000001, 7.000000000000001, 8.0]]
+```
+
 ### Solve a Linear System of Equations
 
 with single right-hand side

Sources/NdArray/LapackError.swift

@@ -4,4 +4,5 @@ public enum LapackError: Error {
  case getrf(_: __CLPK_integer)
  case getri(_: __CLPK_integer)
  case dgesv(_: __CLPK_integer)
+ case gesdd(_: __CLPK_integer)
 }

Sources/NdArray/matrix/Matrix.swift

@@ -513,3 +513,15 @@ public func * (A: Matrix<Float>, B: Matrix<Float>) -> Matrix<Float> {
  cblas_sgemm(order, CblasNoTrans, CblasNoTrans, m, n, k, 1, a.dataStart, lda, b.dataStart, ldb, 0, c.dataStart, ldc)
  return c
 }
+
+public extension Matrix where T: AdditiveArithmetic {
+ convenience init(diag: Vector<T>, order: Contiguous = .C) {
+ let shape = [diag.shape[0], diag.shape[0]]
+ self.init(empty: shape)
+ let z = Vector<T>.zeros(diag.shape[0])
+ for (i, x) in diag.enumerated() {
+ self[Slice(i)] = z
+ self[i, i] = x
+ }
+ }
+}

Sources/NdArray/matrix/svd.swift

@@ -0,0 +1,134 @@
+//
+// Created by Daniel Strobusch on 08.12.22.
+//
+
+import Darwin
+import Accelerate
+
+public extension Matrix where T == Double {
+
+ /// Compute the SVD
+ ///
+ /// Factorization is computed by LAPACKs DGESDD method, which computes a SVD of a general
+ /// M-by-N matrix, optionally computing the left and right singular
+ /// vectors. If singular vectors are desired, it uses a
+ /// divide-and-conquer algorithm.
+ ///
+ /// The SVD is written
+ ///
+ /// A = U * SIGMA * transpose(V)
+ ///
+ /// where SIGMA is an M-by-N matrix which is zero except for its
+ /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
+ /// V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
+ /// are the singular values of A; they are real and non-negative, and
+ /// are returned in descending order. The first min(m,n) columns of
+ /// U and V are the left and right singular vectors of A.
+ ///
+ /// - Returns: a tuple of the matrices (U, SIGMA, transpose(V))
+ func svd() throws -> (Matrix<T>, Vector<T>, Matrix<T>) {
+ // A is destroyed if jobz != O
+ let A = Matrix(copy: self, order: .F)
+
+ var m = __CLPK_integer(shape[0])
+ var n = __CLPK_integer(shape[1])
+
+ // see
+ var jobz: CChar = "A".utf8CString[0]
+
+ // leading dimension is the number of rows in column major order
+ var lda = __CLPK_integer(A.shape[0])
+ var info: __CLPK_integer = 0
+
+ let U: Matrix<T> = Matrix.empty(shape: [Int(m), Int(m)], order: .F)
+ var ldu = __CLPK_integer(U.shape[0])
+ let s: Vector<T> = Vector.empty(shape: [Int(Swift.min(m, n))], order: .F)
+ let Vt: Matrix<T> = Matrix.empty(shape: [Int(n), Int(n)], order: .F)
+ var ldvt = __CLPK_integer(Vt.shape[0])
+
+ let iwork: Vector<__CLPK_integer> = Vector.empty(shape: [8 * Int(Swift.min(m, n))], order: .F)
+
+ // do optimal workspace query
+ var lwork: __CLPK_integer = -1
+ var work: Vector<__CLPK_doublereal> = Vector.empty(shape: [1], order: .F)
+ dgesdd_(&jobz, &m, &n, A.dataStart, &lda, s.dataStart, U.dataStart, &ldu, Vt.dataStart, &ldvt, work.dataStart, &lwork, iwork.dataStart, &info)
+ if info != 0 {
+ throw LapackError.getri(info)
+ }
+
+ // retrieve optimal workspace
+ lwork = __CLPK_integer(work[0])
+ work = Vector.empty(shape: [Int(lwork)], order: .F)
+ dgesdd_(&jobz, &m, &n, A.dataStart, &lda, s.dataStart, U.dataStart, &ldu, Vt.dataStart, &ldvt, work.dataStart, &lwork, iwork.dataStart, &info)
+
+ if info != 0 {
+ throw LapackError.gesdd(info)
+ }
+ return (U, s, Vt)
+ }
+
+}
+
+public extension Matrix where T == Float {
+
+ /// Compute the SVD
+ ///
+ /// Factorization is computed by LAPACKs SGESDD method, which computes a SVD of a general
+ /// M-by-N matrix, optionally computing the left and right singular
+ /// vectors. If singular vectors are desired, it uses a
+ /// divide-and-conquer algorithm.
+ ///
+ /// The SVD is written
+ ///
+ /// A = U * SIGMA * transpose(V)
+ ///
+ /// where SIGMA is an M-by-N matrix which is zero except for its
+ /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
+ /// V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
+ /// are the singular values of A; they are real and non-negative, and
+ /// are returned in descending order. The first min(m,n) columns of
+ /// U and V are the left and right singular vectors of A.
+ ///
+ /// - Returns: a tuple of the matrices (U, SIGMA, transpose(V))
+ func svd() throws -> (Matrix<T>, Vector<T>, Matrix<T>) {
+ // A is destroyed if jobz != O
+ let A = Matrix(copy: self, order: .F)
+
+ var m = __CLPK_integer(shape[0])
+ var n = __CLPK_integer(shape[1])
+
+ // see
+ var jobz: CChar = "A".utf8CString[0]
+
+ // leading dimension is the number of rows in column major order
+ var lda = __CLPK_integer(A.shape[0])
+ var info: __CLPK_integer = 0
+
+ let U: Matrix<T> = Matrix.empty(shape: [Int(m), Int(m)], order: .F)
+ var ldu = __CLPK_integer(U.shape[0])
+ let s: Vector<T> = Vector.empty(shape: [Int(Swift.min(m, n))], order: .F)
+ let Vt: Matrix<T> = Matrix.empty(shape: [Int(n), Int(n)], order: .F)
+ var ldvt = __CLPK_integer(Vt.shape[0])
+
+ let iwork: Vector<__CLPK_integer> = Vector.empty(shape: [8 * Int(Swift.min(m, n))], order: .F)
+
+ // do optimal workspace query
+ var lwork: __CLPK_integer = -1
+ var work: Vector<__CLPK_real> = Vector.empty(shape: [1], order: .F)
+ sgesdd_(&jobz, &m, &n, A.dataStart, &lda, s.dataStart, U.dataStart, &ldu, Vt.dataStart, &ldvt, work.dataStart, &lwork, iwork.dataStart, &info)
+ if info != 0 {
+ throw LapackError.getri(info)
+ }
+
+ // retrieve optimal workspace
+ lwork = __CLPK_integer(work[0])
+ work = Vector.empty(shape: [Int(lwork)], order: .F)
+ sgesdd_(&jobz, &m, &n, A.dataStart, &lda, s.dataStart, U.dataStart, &ldu, Vt.dataStart, &ldvt, work.dataStart, &lwork, iwork.dataStart, &info)
+
+ if info != 0 {
+ throw LapackError.gesdd(info)
+ }
+ return (U, s, Vt)
+ }
+
+}

Tests/NdArrayTests/ReadmeExamples.swift

@@ -178,6 +178,31 @@ class ReadmeExamples: XCTestCase {
  // [[-1.5, 0.5],
  // [ 1.0, 0.0]]
  }
+ do {
+ print("SVD")
+ let A = Matrix<Double>(NdArray.range(from: 1, to: 9).reshaped([2, 4]))
+ let (U, s, Vt) = try A.svd()
+ print(U)
+ // [[-0.3761682344281408, -0.9265513797988838],
+ // [-0.9265513797988838, 0.3761682344281408]]
+ print(s)
+ // [14.227407412633742, 1.2573298353791098]
+ print(Vt)
+ // [[ -0.3520616924890126, -0.44362578258952023, -0.5351898726900277, -0.6267539627905352],
+ // [ 0.7589812676751458, 0.32124159914593237, -0.1164980693832819, -0.554237737912496],
+ // [ -0.4000874340557387, 0.25463292200666415, 0.6909964581538871, -0.5455419461048127],
+ // [ -0.3740722458438949, 0.7969705609558909, -0.471724384380099, 0.04882606926810252]]
+ let Sd = Matrix(diag: s)
+ let S = Matrix<Double>.zeros(A.shape)
+ let mn = A.shape.min()!
+ S[..<mn, ..<mn] = Sd
+ print(S)
+ // [[14.227407412633742, 0.0, 0.0, 0.0],
+ // [ 0.0, 1.2573298353791098, 0.0, 0.0]]
+ print(U * S * Vt)
+ // [[1.0000000000000004, 2.0, 3.0000000000000004, 3.9999999999999996],
+ // [ 4.999999999999999, 6.000000000000001, 7.000000000000001, 8.0]]
+ }
  do {
  print("LU Factorization")
  let A = Matrix<Double>(NdArray.range(to: 4).reshaped([2, 2]))

Tests/NdArrayTests/matrix/svdTestsDouble.swift

@@ -0,0 +1,101 @@
+import XCTest
+import Darwin
+@testable import NdArray
+
+class SvdTestsDouble: XCTestCase {
+ func testSvdUnit() throws {
+ let a = Matrix<Double>([
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ ])
+ let (U, s, Vt) = try a.svd()
+ let S = Matrix(diag: s)
+ XCTAssertEqual((U * S * Vt).dataArray, a.dataArray, accuracy: 1e-15)
+
+ XCTAssertEqual(U.dataArray, Matrix<Double>([
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ ], order: .F).dataArray, accuracy: 1e-15)
+ XCTAssertEqual(S.dataArray, Matrix<Double>([
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ ], order: .F).dataArray, accuracy: 1e-15)
+ XCTAssertEqual(Vt.dataArray, Matrix<Double>([
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ ], order: .F).dataArray, accuracy: 1e-15)
+ }
+
+ func testSvdWide() throws {
+ let a = Matrix<Double>([
+ [1, 2, 3],
+ [4, 5, 6],
+ ])
+ let (U, s, Vt) = try a.svd()
+ let Sd = Matrix(diag: s)
+ let S = Matrix<Double>.zeros(a.shape)
+ let mn = a.shape.min()!
+ S[..<mn, ..<mn] = Sd
+ print("a")
+ print(a)
+ print("U * S * Vt")
+ print(U * S * Vt)
+ XCTAssertEqual(NdArray(U * S * Vt, order: .C).dataArray, a.dataArray, accuracy: 1e-8)
+
+ print("U")
+ print(U)
+ XCTAssertEqual(U.dataArray, Matrix<Double>([
+ [-0.3863177, -0.92236578],
+ [-0.92236578, 0.3863177]
+ ], order: .F).dataArray, accuracy: 1e-8)
+ print("s")
+ print(s)
+ XCTAssertEqual(s.dataArray, Matrix<Double>([[9.508032, 0.77286964]]).dataArray, accuracy: 1e-8)
+ print("Vt")
+ print(Vt)
+ XCTAssertEqual(Vt.dataArray, Matrix<Double>([
+ [-0.42866713, -0.56630692, -0.7039467],
+ [0.80596391, 0.11238241, -0.58119908],
+ [0.40824829, -0.81649658, 0.40824829],
+ ], order: .F).dataArray, accuracy: 1e-8)
+ }
+
+ func testSvdTall() throws {
+ let a = Matrix<Double>([
+ [1, 2],
+ [3, 4],
+ [5, 6],
+ ])
+ let (U, s, Vt) = try a.svd()
+ let Sd = Matrix(diag: s)
+ let S = Matrix<Double>.zeros(a.shape)
+ let mn = a.shape.min()!
+ S[..<mn, ..<mn] = Sd
+ print("a")
+ print(a)
+ print("U * S * Vt")
+ print(U * S * Vt)
+ XCTAssertEqual(NdArray(U * S * Vt, order: .C).dataArray, a.dataArray, accuracy: 1e-8)
+
+ print("U")
+ print(U)
+ XCTAssertEqual(U.dataArray, Matrix<Double>([
+ [-0.2298477, 0.88346102, 0.40824829],
+ [-0.52474482, 0.24078249, -0.81649658],
+ [-0.81964194, -0.40189603, 0.40824829]
+ ], order: .F).dataArray, accuracy: 1e-8)
+ print("s")
+ print(s)
+ XCTAssertEqual(s.dataArray, Matrix<Double>([[9.52551809, 0.51430058]]).dataArray, accuracy: 1e-8)
+ print("Vt")
+ print(Vt)
+ XCTAssertEqual(Vt.dataArray, Matrix<Double>([
+ [-0.61962948, -0.78489445],
+ [-0.78489445, 0.61962948],
+ ], order: .F).dataArray, accuracy: 1e-8)
+ }
+}