diff --git a/README.md b/README.md index 359240c..a06db8c 100644 --- a/README.md +++ b/README.md @@ -7,8 +7,8 @@ by Accelerate. This package should provide small, explicit Swift APIs for selecting high-volume row submatrices from tall dense matrices. The first implementation target is a -real-valued `Double` path using Accelerate BLAS/LAPACK, with `Float` support -following the same shape once the API is settled. +real-valued `Double` and `Float` path using Accelerate BLAS/LAPACK without +deprecated compatibility shims. ## Implementation Plan @@ -17,11 +17,11 @@ following the same shape once the API is settled. - Keep dimensions explicit and validate shape before calling Accelerate. - Return descriptive errors for non-tall, rank-deficient, or malformed input. -2. Implement square MaxVol for `Double`. +2. Implement square MaxVol for `Double` and `Float`. - Accept an `N x r` matrix where `N >= r`. - - Initialize pivots with `dgetrf`. + - Initialize pivots with `dgetrf` / `sgetrf`. - Solve for expansion coefficients with triangular solves. - - Iterate row swaps using the standard rank-one update with `dger`. + - Iterate row swaps using the standard rank-one update with `dger` / `sger`. - Return selected row indices, coefficient matrix, iteration count, and whether the coefficient tolerance was reached. @@ -41,7 +41,8 @@ following the same shape once the API is settled. 5. Add performance-focused refinements only after the behavior is stable. - Avoid repeated temporary allocations in the swap loop. - Reuse workspace buffers. - - Add `Float` variants mapped to `sgetrf`, `strsm`, and `sger`. + - Keep `Double` and `Float` paths mapped to the matching Accelerate entry + points. - Consider complex support only after real-valued APIs are stable. ## Non-Goals For The First Pass diff --git a/ROADMAP.md b/ROADMAP.md index a4a3a72..069c236 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -2,17 +2,19 @@ ## Current Focus -- Build from `v0.5.0` toward `v1.0.0` with real-valued `Double` algorithms - backed by modern Accelerate BLAS/LAPACK. +- Build from `v0.7.0` toward `v1.0.0` with real-valued `Double` and `Float` + algorithms backed by modern Accelerate BLAS/LAPACK. - Keep Swift Testing as the package test surface and require reference fixtures before broadening the algorithm surface. -- Treat public API compatibility as pre-1.0 until RectMaxVol, Float support, and - reference parity are complete. +- Treat public API compatibility as pre-1.0 until reference parity, + performance behavior, and documentation are complete enough to support a + stable contract. ## Algorithm Milestones -1. Implement square MaxVol for `Double`. - - Accept a tall or square `DenseColumnMajorMatrix` where +1. Implement square MaxVol for `Double` and `Float`. + - Accept a tall or square `DenseColumnMajorMatrix` or + `DenseColumnMajorMatrix` where `rows >= columns`. - Initialize pivots with Accelerate LAPACK LU factorization. - Solve expansion coefficients with Accelerate triangular solve routines. @@ -20,24 +22,24 @@ - Return selected row indices, coefficient matrix, iteration count, and convergence status. -2. Add RectMaxVol for `Double`. +2. Add RectMaxVol for `Double` and `Float`. - Start from square MaxVol selected rows. - Append rows while coefficient row norms exceed tolerance. - Support explicit `minRows` and `maxRows` bounds. - Reuse coefficient-update primitives from square MaxVol where practical. -3. Add `Float` support after the `Double` API is stable. - - Mirror the `Double` API shape. - - Map to the matching single-precision Accelerate entry points. - - Keep shared validation and result semantics identical. - -4. Add reference-parity fixtures for each supported algorithm. +3. Add reference-parity fixtures for each supported algorithm. - Keep small deterministic fixtures generated from `maxvolpy`, `Maxvol.jl`, or an R implementation. - Record selected rows, coefficients, iteration counts, and convergence status for fixture matrices that require zero, one, and multiple swaps. - Include independent reconstruction checks for every fixture. +4. Add performance-focused refinements. + - Reduce temporary allocations in square row swaps and rectangular appends. + - Reuse workspace buffers where it improves measured throughput. + - Add benchmark coverage for representative row counts and ranks. + ## Test Coverage - Cover matrix storage, checked access, shape validation, and result invariants. @@ -47,8 +49,9 @@ maximum-iteration cases. - Keep reference fixtures generated from upstream Python, Julia, or R implementations when algorithm behavior changes. -- Add randomized orthonormal-matrix tests similar to `Maxvol.jl` before `1.0.0`. -- Add Release-mode validation once behavior depends on optimization-sensitive +- Keep randomized orthonormal-matrix tests similar to `Maxvol.jl` across + supported scalar types. +- Keep Release-mode validation in the release path for optimization-sensitive Accelerate calls. - Keep validation clean under the modern Accelerate `ACCELERATE_NEW_LAPACK` and `ACCELERATE_LAPACK_ILP64` import surface. @@ -61,8 +64,8 @@ - Document square MaxVol and RectMaxVol usage with small examples. - Keep algorithm limitations and non-goals visible before `1.0.0`. - Keep public API docs aligned with tested behavior. -- Add a DocC article that explains tolerance, convergence status, and - iteration-limited partial results before `1.0.0`. +- Expand DocC with algorithm notes, limitations, and reference-fixture + provenance before `1.0.0`. ## Swift Package Index @@ -81,14 +84,14 @@ ## Before `1.0.0` -- Complete RectMaxVol for `Double` with reference-parity fixtures. -- Add `Float` support with the same API shape and reference fixtures. - Decide whether complex-valued matrices are in scope for `1.0.0` or explicitly post-1.0. -- Add broader randomized numerical tests, including orthonormal tall matrices - and near-rank-deficient cases. -- Add Release-mode validation for optimization-sensitive Accelerate behavior. -- Add performance benchmarks for allocation count and row-swap throughput. -- Expand DocC with algorithm notes, limitations, and reference-fixture - provenance. +- Add broader randomized numerical tests, including stress cases for larger + ranks, wider condition-number ranges, and deterministic tie behavior. +- Add performance benchmarks for allocation count, row-swap throughput, and + rectangular append throughput. +- Profile whether reusable workspaces are worth adding to the public API or + should remain an internal optimization. +- Expand DocC with reference-fixture provenance and a short comparison with + Python, Julia, and R implementations. - Verify Swift Package Index renders the tagged `v1.0.0` documentation cleanly. diff --git a/Sources/MaxVol/MaxVol.docc/MaxVol.md b/Sources/MaxVol/MaxVol.docc/MaxVol.md index 672d47a..6761a46 100644 --- a/Sources/MaxVol/MaxVol.docc/MaxVol.md +++ b/Sources/MaxVol/MaxVol.docc/MaxVol.md @@ -6,7 +6,7 @@ Accelerate. ## Overview MaxVol provides Swift-native APIs for selecting representative rows from tall -dense matrices. The first implementation path supports real-valued `Double` +dense matrices. The implementation supports real-valued `Double` and `Float` matrices stored in column-major order so the package can call Accelerate BLAS and LAPACK routines directly. @@ -26,7 +26,12 @@ let matrix = try DenseColumnMajorMatrix( ) ``` -Call ``maxVol(_:options:)`` to select a square basis from a tall matrix: +Use `DenseColumnMajorMatrix` for the single-precision path; the public +algorithm calls keep the same shape. + +Call ``maxVol(_:options:)->MaxVolResult`` or +``maxVol(_:options:)->MaxVolResult`` to select a square basis from a +tall matrix: ```swift let result = try maxVol(matrix) @@ -44,8 +49,9 @@ rows: let coefficients = result.coefficients ``` -Call ``rectMaxVol(_:options:)`` when the basis may contain more rows than the -matrix column count: +Call ``rectMaxVol(_:options:)->MaxVolResult`` or +``rectMaxVol(_:options:)->MaxVolResult`` when the basis may contain more +rows than the matrix column count: ```swift let rectangular = try rectMaxVol(matrix, options: RectMaxVolOptions(minRows: 3)) @@ -59,8 +65,10 @@ let rectangular = try rectMaxVol(matrix, options: RectMaxVolOptions(minRows: 3)) ### MaxVol -- ``maxVol(_:options:)`` -- ``rectMaxVol(_:options:)`` +- ``maxVol(_:options:)->MaxVolResult`` +- ``maxVol(_:options:)->MaxVolResult`` +- ``rectMaxVol(_:options:)->MaxVolResult`` +- ``rectMaxVol(_:options:)->MaxVolResult`` - ``MaxVolOptions`` - ``RectMaxVolOptions`` - ``MaxVolResult`` diff --git a/Sources/MaxVol/MaxVol.docc/ToleranceAndConvergence.md b/Sources/MaxVol/MaxVol.docc/ToleranceAndConvergence.md index 6f2a527..5f59734 100644 --- a/Sources/MaxVol/MaxVol.docc/ToleranceAndConvergence.md +++ b/Sources/MaxVol/MaxVol.docc/ToleranceAndConvergence.md @@ -22,9 +22,10 @@ let coefficients = result.coefficients ## Square MaxVol -``maxVol(_:options:)`` selects exactly one row per matrix column. The square -algorithm starts from an LU-pivoted basis, computes expansion coefficients, and -swaps rows while any coefficient magnitude is larger than +``maxVol(_:options:)->MaxVolResult`` and +``maxVol(_:options:)->MaxVolResult`` select exactly one row per matrix +column. The square algorithm starts from an LU-pivoted basis, computes expansion +coefficients, and swaps rows while any coefficient magnitude is larger than ``MaxVolOptions/tolerance``. The square tolerance must be at least `1.0`. Values closer to `1.0` usually do @@ -37,10 +38,11 @@ validated and reconstructs through its current selected rows, but ## RectMaxVol -``rectMaxVol(_:options:)`` starts from square MaxVol, then appends extra rows. -The rectangular stopping test uses coefficient row norms instead of individual -coefficient magnitudes: unselected rows are appended while their coefficient row -norm exceeds ``RectMaxVolOptions/tolerance``. +``rectMaxVol(_:options:)->MaxVolResult`` and +``rectMaxVol(_:options:)->MaxVolResult`` start from square MaxVol, then +append extra rows. The rectangular stopping test uses coefficient row norms +instead of individual coefficient magnitudes: unselected rows are appended while +their coefficient row norm exceeds ``RectMaxVolOptions/tolerance``. ``RectMaxVolOptions/minRows`` can force extra rows even when the tolerance is already satisfied. ``RectMaxVolOptions/maxRows`` can stop the append loop before @@ -57,12 +59,18 @@ stronger square basis matters. For rectangular selection, set basis size, and set ``RectMaxVolOptions/maxRows`` when runtime or storage must be bounded. +The same option types apply to both `Double` and `Float` overloads. Tolerances +are expressed as `Double` values so callers can keep one configuration surface +while the numerical work uses the matrix scalar type. + ## Topics ### Related APIs -- ``maxVol(_:options:)`` -- ``rectMaxVol(_:options:)`` +- ``maxVol(_:options:)->MaxVolResult`` +- ``maxVol(_:options:)->MaxVolResult`` +- ``rectMaxVol(_:options:)->MaxVolResult`` +- ``rectMaxVol(_:options:)->MaxVolResult`` - ``MaxVolOptions`` - ``RectMaxVolOptions`` - ``MaxVolResult`` diff --git a/Sources/MaxVol/MaxVol.swift b/Sources/MaxVol/MaxVol.swift index 065e899..06ff41a 100644 --- a/Sources/MaxVol/MaxVol.swift +++ b/Sources/MaxVol/MaxVol.swift @@ -1,7 +1,3 @@ -import Accelerate - -private typealias LAPACKInt = __LAPACK_int - /// Selects a high-volume square row basis from a tall dense `Double` matrix. /// /// The returned coefficients are shaped so the input matrix `A` can be @@ -11,6 +7,25 @@ public func maxVol( _ matrix: DenseColumnMajorMatrix, options: MaxVolOptions = MaxVolOptions() ) throws -> MaxVolResult { + try maxVolImpl(matrix, options: options) +} + +/// Selects a high-volume square row basis from a tall dense `Float` matrix. +/// +/// The returned coefficients are shaped so the input matrix `A` can be +/// approximated by `C * A[selectedRows, :]`, where `C` is +/// ``MaxVolResult/coefficients``. +public func maxVol( + _ matrix: DenseColumnMajorMatrix, + options: MaxVolOptions = MaxVolOptions() +) throws -> MaxVolResult { + try maxVolImpl(matrix, options: options) +} + +func maxVolImpl( + _ matrix: DenseColumnMajorMatrix, + options: MaxVolOptions = MaxVolOptions() +) throws -> MaxVolResult { let input = try matrix.validatedTallMatrix() let options = try options.validated() var selectedRows = try initialPivotRows(for: input) @@ -52,7 +67,9 @@ private struct CoefficientPivot { let value: Double } -private func initialPivotRows(for matrix: DenseColumnMajorMatrix) throws -> [Int] { +private func initialPivotRows( + for matrix: DenseColumnMajorMatrix +) throws -> [Int] { var factorization = matrix.values var rowCount = try lapackInt(matrix.rows) var columnCount = try lapackInt(matrix.columns) @@ -62,18 +79,18 @@ private func initialPivotRows(for matrix: DenseColumnMajorMatrix) throws factorization.withUnsafeMutableBufferPointer { factorBuffer -> Void in pivots.withUnsafeMutableBufferPointer { pivotBuffer -> Void in - dgetrf_( - &rowCount, - &columnCount, - factorBuffer.baseAddress, - &leadingDimension, - pivotBuffer.baseAddress, - &info + Scalar.getrf( + rowCount: &rowCount, + columnCount: &columnCount, + values: factorBuffer.baseAddress, + leadingDimension: &leadingDimension, + pivots: pivotBuffer.baseAddress, + info: &info ) } } - try validateLapackInfo(info, routine: "dgetrf", rankDeficientInfoIsPivot: true) + try validateLapackInfo(info, routine: Scalar.getrfRoutineName, rankDeficientInfoIsPivot: true) var permutation = Array(0..) throws return Array(permutation.prefix(matrix.columns)) } -private func expansionCoefficients( - for matrix: DenseColumnMajorMatrix, +private func expansionCoefficients( + for matrix: DenseColumnMajorMatrix, selectedRows: [Int] -) throws -> DenseColumnMajorMatrix { +) throws -> DenseColumnMajorMatrix { let rank = selectedRows.count var basisValues = selectedRowsForBasis(matrix: matrix, selectedRows: selectedRows) var rowRankDimension = try lapackInt(rank) @@ -98,18 +115,18 @@ private func expansionCoefficients( basisValues.withUnsafeMutableBufferPointer { basisBuffer -> Void in pivots.withUnsafeMutableBufferPointer { pivotBuffer -> Void in - dgetrf_( - &rowRankDimension, - &columnRankDimension, - basisBuffer.baseAddress, - &leadingDimension, - pivotBuffer.baseAddress, - &info + Scalar.getrf( + rowCount: &rowRankDimension, + columnCount: &columnRankDimension, + values: basisBuffer.baseAddress, + leadingDimension: &leadingDimension, + pivots: pivotBuffer.baseAddress, + info: &info ) } } - try validateLapackInfo(info, routine: "dgetrf", rankDeficientInfoIsPivot: true) + try validateLapackInfo(info, routine: Scalar.getrfRoutineName, rankDeficientInfoIsPivot: true) try validateNonsingularFactorization(basisValues, dimension: rank) var transposedRightHandSide = transposedValues(matrix) @@ -123,22 +140,22 @@ private func expansionCoefficients( basisValues.withUnsafeMutableBufferPointer { basisBuffer -> Void in pivots.withUnsafeMutableBufferPointer { pivotBuffer -> Void in transposedRightHandSide.withUnsafeMutableBufferPointer { rightHandSideBuffer -> Void in - dgetrs_( - &transpose, - &solveRankDimension, - &rightHandSides, - basisBuffer.baseAddress, - &basisLeadingDimension, - pivotBuffer.baseAddress, - rightHandSideBuffer.baseAddress, - &rightHandSideLeadingDimension, - &info + Scalar.getrs( + transpose: &transpose, + dimension: &solveRankDimension, + rightHandSides: &rightHandSides, + factorization: basisBuffer.baseAddress, + leadingDimension: &basisLeadingDimension, + pivots: pivotBuffer.baseAddress, + rightHandSide: rightHandSideBuffer.baseAddress, + rightHandSideLeadingDimension: &rightHandSideLeadingDimension, + info: &info ) } } } - try validateLapackInfo(info, routine: "dgetrs", rankDeficientInfoIsPivot: false) + try validateLapackInfo(info, routine: Scalar.getrsRoutineName, rankDeficientInfoIsPivot: false) let coefficientValues = (0.., +private func selectedRowsForBasis( + matrix: DenseColumnMajorMatrix, selectedRows: [Int] -) -> [Double] { +) -> [Scalar] { (0..) -> [Double] { +private func transposedValues( + _ matrix: DenseColumnMajorMatrix +) -> [Scalar] { (0..) -> [Doub } } -private func maximumMagnitude(in coefficients: DenseColumnMajorMatrix) -> CoefficientPivot { - var pivot = CoefficientPivot(row: 0, column: 0, value: abs(coefficients[row: 0, column: 0])) +private func maximumMagnitude( + in coefficients: DenseColumnMajorMatrix +) -> CoefficientPivot { + var pivot = CoefficientPivot( + row: 0, + column: 0, + value: coefficients[row: 0, column: 0].magnitudeAsDouble + ) for column in 0.. pivot.value { pivot = CoefficientPivot(row: row, column: column, value: magnitude) } @@ -187,10 +212,10 @@ private func maximumMagnitude(in coefficients: DenseColumnMajorMatrix) - return pivot } -private func replaceBasisRow( +private func replaceBasisRow( pivotRow: Int, pivotColumn: Int, - coefficients: inout DenseColumnMajorMatrix + coefficients: inout DenseColumnMajorMatrix ) throws { let gamma = coefficients[row: pivotRow, column: pivotColumn] guard gamma != 0 else { @@ -209,17 +234,16 @@ private func replaceBasisRow( replacementColumn.withUnsafeBufferPointer { columnBuffer -> Void in replacementRow.withUnsafeBufferPointer { rowBuffer -> Void in coefficients.values.withUnsafeMutableBufferPointer { coefficientBuffer -> Void in - cblas_dger( - CblasColMajor, - rowCount, - columnCount, - -1, - columnBuffer.baseAddress, - increment, - rowBuffer.baseAddress, - increment, - coefficientBuffer.baseAddress, - leadingDimension + Scalar.rankOneUpdate( + rowCount: rowCount, + columnCount: columnCount, + alpha: -1, + x: columnBuffer.baseAddress, + incrementX: increment, + y: rowBuffer.baseAddress, + incrementY: increment, + values: coefficientBuffer.baseAddress, + leadingDimension: leadingDimension ) } } @@ -230,38 +254,17 @@ private func replaceBasisRow( } } -private func lapackInt(_ value: Int) throws -> LAPACKInt { - guard value <= Int(LAPACKInt.max) else { - throw MaxVolError.invalidDimensions(rows: value, columns: value) - } - - return LAPACKInt(value) -} - -private func validateNonsingularFactorization(_ values: [Double], dimension: Int) throws { - let scale = max(values.map(abs).max() ?? 0, 1) - let threshold = Double.ulpOfOne * Double(dimension) * scale +private func validateNonsingularFactorization( + _ values: [Scalar], + dimension: Int +) throws { + let scale = max(values.map(\.magnitudeAsDouble).max() ?? 0, 1) + let threshold = Scalar.rankToleranceUnit * Double(dimension) * scale for pivot in 0.. threshold else { + guard diagonal.magnitudeAsDouble > threshold else { throw MaxVolError.rankDeficient(pivot: pivot + 1) } } } - -private func validateLapackInfo( - _ info: LAPACKInt, - routine: String, - rankDeficientInfoIsPivot: Bool -) throws { - if info < 0 { - throw MaxVolError.lapackFailure(routine: routine, info: Int(info)) - } - if info > 0 { - if rankDeficientInfoIsPivot { - throw MaxVolError.rankDeficient(pivot: Int(info)) - } - throw MaxVolError.lapackFailure(routine: routine, info: Int(info)) - } -} diff --git a/Sources/MaxVol/MaxVolScalar.swift b/Sources/MaxVol/MaxVolScalar.swift new file mode 100644 index 0000000..1421927 --- /dev/null +++ b/Sources/MaxVol/MaxVolScalar.swift @@ -0,0 +1,218 @@ +import Accelerate + +typealias LAPACKInt = __LAPACK_int + +protocol MaxVolScalar: BinaryFloatingPoint, Sendable { + static var getrfRoutineName: String { get } + static var getrsRoutineName: String { get } + static var rankToleranceUnit: Double { get } + + var magnitudeAsDouble: Double { get } + + static func getrf( + rowCount: inout LAPACKInt, + columnCount: inout LAPACKInt, + values: UnsafeMutablePointer?, + leadingDimension: inout LAPACKInt, + pivots: UnsafeMutablePointer?, + info: inout LAPACKInt + ) + + static func getrs( + transpose: inout CChar, + dimension: inout LAPACKInt, + rightHandSides: inout LAPACKInt, + factorization: UnsafeMutablePointer?, + leadingDimension: inout LAPACKInt, + pivots: UnsafeMutablePointer?, + rightHandSide: UnsafeMutablePointer?, + rightHandSideLeadingDimension: inout LAPACKInt, + info: inout LAPACKInt + ) + + static func rankOneUpdate( + rowCount: LAPACKInt, + columnCount: LAPACKInt, + alpha: Self, + x: UnsafePointer?, + incrementX: LAPACKInt, + y: UnsafePointer?, + incrementY: LAPACKInt, + values: UnsafeMutablePointer?, + leadingDimension: LAPACKInt + ) +} + +extension Double: MaxVolScalar { + static var getrfRoutineName: String { "dgetrf" } + static var getrsRoutineName: String { "dgetrs" } + static var rankToleranceUnit: Double { Double.ulpOfOne } + + var magnitudeAsDouble: Double { abs(self) } + + static func getrf( + rowCount: inout LAPACKInt, + columnCount: inout LAPACKInt, + values: UnsafeMutablePointer?, + leadingDimension: inout LAPACKInt, + pivots: UnsafeMutablePointer?, + info: inout LAPACKInt + ) { + dgetrf_( + &rowCount, + &columnCount, + values, + &leadingDimension, + pivots, + &info + ) + } + + static func getrs( + transpose: inout CChar, + dimension: inout LAPACKInt, + rightHandSides: inout LAPACKInt, + factorization: UnsafeMutablePointer?, + leadingDimension: inout LAPACKInt, + pivots: UnsafeMutablePointer?, + rightHandSide: UnsafeMutablePointer?, + rightHandSideLeadingDimension: inout LAPACKInt, + info: inout LAPACKInt + ) { + dgetrs_( + &transpose, + &dimension, + &rightHandSides, + factorization, + &leadingDimension, + pivots, + rightHandSide, + &rightHandSideLeadingDimension, + &info + ) + } + + static func rankOneUpdate( + rowCount: LAPACKInt, + columnCount: LAPACKInt, + alpha: Double, + x: UnsafePointer?, + incrementX: LAPACKInt, + y: UnsafePointer?, + incrementY: LAPACKInt, + values: UnsafeMutablePointer?, + leadingDimension: LAPACKInt + ) { + cblas_dger( + CblasColMajor, + rowCount, + columnCount, + alpha, + x, + incrementX, + y, + incrementY, + values, + leadingDimension + ) + } +} + +extension Float: MaxVolScalar { + static var getrfRoutineName: String { "sgetrf" } + static var getrsRoutineName: String { "sgetrs" } + static var rankToleranceUnit: Double { Double(Float.ulpOfOne) } + + var magnitudeAsDouble: Double { Double(abs(self)) } + + static func getrf( + rowCount: inout LAPACKInt, + columnCount: inout LAPACKInt, + values: UnsafeMutablePointer?, + leadingDimension: inout LAPACKInt, + pivots: UnsafeMutablePointer?, + info: inout LAPACKInt + ) { + sgetrf_( + &rowCount, + &columnCount, + values, + &leadingDimension, + pivots, + &info + ) + } + + static func getrs( + transpose: inout CChar, + dimension: inout LAPACKInt, + rightHandSides: inout LAPACKInt, + factorization: UnsafeMutablePointer?, + leadingDimension: inout LAPACKInt, + pivots: UnsafeMutablePointer?, + rightHandSide: UnsafeMutablePointer?, + rightHandSideLeadingDimension: inout LAPACKInt, + info: inout LAPACKInt + ) { + sgetrs_( + &transpose, + &dimension, + &rightHandSides, + factorization, + &leadingDimension, + pivots, + rightHandSide, + &rightHandSideLeadingDimension, + &info + ) + } + + static func rankOneUpdate( + rowCount: LAPACKInt, + columnCount: LAPACKInt, + alpha: Float, + x: UnsafePointer?, + incrementX: LAPACKInt, + y: UnsafePointer?, + incrementY: LAPACKInt, + values: UnsafeMutablePointer?, + leadingDimension: LAPACKInt + ) { + cblas_sger( + CblasColMajor, + rowCount, + columnCount, + alpha, + x, + incrementX, + y, + incrementY, + values, + leadingDimension + ) + } +} + +func lapackInt(_ value: Int) throws -> LAPACKInt { + guard value <= Int(LAPACKInt.max) else { + throw MaxVolError.invalidDimensions(rows: value, columns: value) + } + + return LAPACKInt(value) +} + +func validateLapackInfo( + _ info: LAPACKInt, + routine: String, + rankDeficientInfoIsPivot: Bool +) throws { + if info < 0 { + throw MaxVolError.lapackFailure(routine: routine, info: Int(info)) + } + if info > 0 { + if rankDeficientInfoIsPivot { + throw MaxVolError.rankDeficient(pivot: Int(info)) + } + throw MaxVolError.lapackFailure(routine: routine, info: Int(info)) + } +} diff --git a/Sources/MaxVol/RectMaxVol.swift b/Sources/MaxVol/RectMaxVol.swift index 2c233f5..f88450b 100644 --- a/Sources/MaxVol/RectMaxVol.swift +++ b/Sources/MaxVol/RectMaxVol.swift @@ -1,19 +1,34 @@ -import Accelerate - -private typealias LAPACKInt = __LAPACK_int - /// Selects a high-volume rectangular row basis from a tall dense `Double` matrix. /// -/// RectMaxVol starts with ``maxVol(_:options:)`` and appends rows until every -/// remaining unselected coefficient row satisfies ``RectMaxVolOptions/tolerance`` -/// or the configured row bounds stop the append loop. +/// RectMaxVol starts with ``maxVol(_:options:)->MaxVolResult`` and appends +/// rows until every remaining unselected coefficient row satisfies +/// ``RectMaxVolOptions/tolerance`` or the configured row bounds stop the append loop. public func rectMaxVol( _ matrix: DenseColumnMajorMatrix, options: RectMaxVolOptions = RectMaxVolOptions() ) throws -> MaxVolResult { + try rectMaxVolImpl(matrix, options: options) +} + +/// Selects a high-volume rectangular row basis from a tall dense `Float` matrix. +/// +/// RectMaxVol starts with ``maxVol(_:options:)->MaxVolResult`` and appends +/// rows until every remaining unselected coefficient row satisfies +/// ``RectMaxVolOptions/tolerance`` or the configured row bounds stop the append loop. +public func rectMaxVol( + _ matrix: DenseColumnMajorMatrix, + options: RectMaxVolOptions = RectMaxVolOptions() +) throws -> MaxVolResult { + try rectMaxVolImpl(matrix, options: options) +} + +func rectMaxVolImpl( + _ matrix: DenseColumnMajorMatrix, + options: RectMaxVolOptions = RectMaxVolOptions() +) throws -> MaxVolResult { let input = try matrix.validatedTallMatrix() - let options = try options.resolved(for: input) - let initial = try maxVol( + let options = try options.resolved(rows: input.rows, columns: input.columns) + let initial = try maxVolImpl( input, options: MaxVolOptions(maxIterations: options.startMaxVolIterations) ) @@ -69,8 +84,8 @@ private struct RowNormCandidate { let value: Double } -private func maximumUnselectedRowNormSquared( - in coefficients: DenseColumnMajorMatrix, +private func maximumUnselectedRowNormSquared( + in coefficients: DenseColumnMajorMatrix, selectedRows: [Int] ) -> RowNormCandidate { let selected = Set(selectedRows) @@ -79,7 +94,7 @@ private func maximumUnselectedRowNormSquared( for row in 0.. candidate.value { @@ -90,13 +105,13 @@ private func maximumUnselectedRowNormSquared( return candidate } -private func appendRectangularBasisRow( +private func appendRectangularBasisRow( _ candidateRow: Int, - to coefficients: DenseColumnMajorMatrix -) throws -> DenseColumnMajorMatrix { + to coefficients: DenseColumnMajorMatrix +) throws -> DenseColumnMajorMatrix { let candidateCoefficients = try coefficients.row(candidateRow) let projection = (0.. Void in candidateCoefficients.withUnsafeBufferPointer { coefficientBuffer -> Void in updatedValues.withUnsafeMutableBufferPointer { updatedBuffer -> Void in - cblas_dger( - CblasColMajor, - rowCount, - columnCount, - -scale, - projectionBuffer.baseAddress, - increment, - coefficientBuffer.baseAddress, - increment, - updatedBuffer.baseAddress, - leadingDimension + Scalar.rankOneUpdate( + rowCount: rowCount, + columnCount: columnCount, + alpha: -scale, + x: projectionBuffer.baseAddress, + incrementX: increment, + y: coefficientBuffer.baseAddress, + incrementY: increment, + values: updatedBuffer.baseAddress, + leadingDimension: leadingDimension ) } } @@ -135,10 +149,10 @@ private func appendRectangularBasisRow( ) } -private func coefficientsWithIdentityRows( - _ coefficients: DenseColumnMajorMatrix, +private func coefficientsWithIdentityRows( + _ coefficients: DenseColumnMajorMatrix, selectedRows: [Int] -) throws -> DenseColumnMajorMatrix { +) throws -> DenseColumnMajorMatrix { var output = coefficients for (identityColumn, selectedRow) in selectedRows.enumerated() { @@ -153,11 +167,3 @@ private func coefficientsWithIdentityRows( return output } - -private func lapackInt(_ value: Int) throws -> LAPACKInt { - guard value <= Int(LAPACKInt.max) else { - throw MaxVolError.invalidDimensions(rows: value, columns: value) - } - - return LAPACKInt(value) -} diff --git a/Sources/MaxVol/RectMaxVolOptions.swift b/Sources/MaxVol/RectMaxVolOptions.swift index 729b855..5be040b 100644 --- a/Sources/MaxVol/RectMaxVolOptions.swift +++ b/Sources/MaxVol/RectMaxVolOptions.swift @@ -32,7 +32,11 @@ public struct RectMaxVolOptions: Equatable, Hashable, Sendable { self.startMaxVolIterations = startMaxVolIterations } - func resolved(for matrix: DenseColumnMajorMatrix) throws -> ResolvedRectMaxVolOptions { + func resolved(for matrix: DenseColumnMajorMatrix) throws -> ResolvedRectMaxVolOptions { + try resolved(rows: matrix.rows, columns: matrix.columns) + } + + func resolved(rows: Int, columns: Int) throws -> ResolvedRectMaxVolOptions { guard tolerance.isFinite, tolerance > 0 else { throw MaxVolError.invalidTolerance(tolerance) } @@ -40,21 +44,21 @@ public struct RectMaxVolOptions: Equatable, Hashable, Sendable { throw MaxVolError.invalidIterationLimit(startMaxVolIterations) } - let requiredRows = matrix.columns + let requiredRows = columns let resolvedMinRows = minRows ?? requiredRows - let resolvedMaxRows = maxRows ?? matrix.rows + let resolvedMaxRows = maxRows ?? rows guard resolvedMinRows >= requiredRows, resolvedMaxRows >= requiredRows, resolvedMinRows <= resolvedMaxRows, - resolvedMaxRows <= matrix.rows + resolvedMaxRows <= rows else { throw MaxVolError.invalidRowSelectionBounds( minRows: resolvedMinRows, maxRows: resolvedMaxRows, requiredRows: requiredRows, - availableRows: matrix.rows + availableRows: rows ) } diff --git a/Tests/MaxVolTests/FloatMaxVolTests.swift b/Tests/MaxVolTests/FloatMaxVolTests.swift new file mode 100644 index 0000000..71e37c3 --- /dev/null +++ b/Tests/MaxVolTests/FloatMaxVolTests.swift @@ -0,0 +1,391 @@ +@testable import MaxVol +import Testing + +struct FloatMaxVolTests { + // Reference fixture values in this suite were generated with + // maxvolpy.maxvol.py_maxvol and py_rect_maxvol running on NumPy/SciPy + // with np.float32 inputs. + + @Test func squareReferenceFixtureMatchesFloat32PivotsAndCoefficients() throws { + let matrix = try DenseColumnMajorMatrix( + rows: 3, + columns: 2, + rowMajorValues: [ + 0.25, 0.0, + 0.5, 1.0, + 0.5, 1.5, + ] + ) + + let result = try maxVol(matrix, options: MaxVolOptions(tolerance: 1.0)) + + #expect(result.selectedRows == [2, 0]) + #expect(result.iterations == 1) + #expect(result.converged) + try expectCoefficients( + result.coefficients, + rowMajorValues: [ + 0.0, 1.0, + 0.666_666_626_930_236_8, 0.666_666_686_534_881_6, + 1.0, 0.0, + ] + ) + #expect(maximumAbsoluteCoefficient(in: result.coefficients) == 1.0) + try expectReconstruction(of: matrix, using: result) + } + + @Test func zeroMaxIterationsReturnsFloat32InitialLUFixture() throws { + let matrix = try DenseColumnMajorMatrix( + rows: 3, + columns: 2, + rowMajorValues: [ + 0.25, 0.0, + 0.5, 1.0, + 0.5, 1.5, + ] + ) + + let result = try maxVol(matrix, options: MaxVolOptions(maxIterations: 0)) + + #expect(result.selectedRows == [1, 0]) + #expect(result.iterations == 0) + #expect(!result.converged) + try expectCoefficients( + result.coefficients, + rowMajorValues: [ + 0.0, 1.0, + 1.0, 0.0, + 1.5, -1.0, + ] + ) + #expect(maximumAbsoluteCoefficient(in: result.coefficients) == 1.5) + try expectReconstruction(of: matrix, using: result) + } + + @Test func iterationLimitReturnsFloat32PartialFixture() throws { + let matrix = try DenseColumnMajorMatrix( + rows: 4, + columns: 2, + rowMajorValues: [ + 0.25, 0.0, + 0.5, 1.0, + 0.5, 0.5, + 0.5, 1.5, + ] + ) + + let limited = try maxVol(matrix, options: MaxVolOptions(maxIterations: 1)) + let converged = try maxVol(matrix, options: MaxVolOptions(maxIterations: 2)) + + #expect(limited.selectedRows == [3, 0]) + #expect(limited.iterations == 1) + #expect(!limited.converged) + try expectCoefficients( + limited.coefficients, + rowMajorValues: [ + 0.0, 1.0, + 0.666_666_626_930_236_8, 0.666_666_686_534_881_6, + 0.333_333_313_465_118_4, 1.333_333_373_069_763_2, + 1.0, 0.0, + ] + ) + try expectReconstruction(of: matrix, using: limited) + + #expect(converged.selectedRows == [3, 2]) + #expect(converged.iterations == 2) + #expect(converged.converged) + try expectCoefficients( + converged.coefficients, + rowMajorValues: [ + -0.249_999_985_098_838_8, 0.75, + 0.499_999_970_197_677_6, 0.5, + 0.0, 1.0, + 1.0, 0.0, + ] + ) + try expectReconstruction(of: matrix, using: converged) + } + + @Test func rectangularReferenceFixtureAppendsRequiredMinimumRows() throws { + let matrix = try DenseColumnMajorMatrix( + rows: 4, + columns: 2, + rowMajorValues: [ + 1.0, 0.0, + 0.0, 1.0, + 0.5, 0.25, + -0.25, 0.75, + ] + ) + + let result = try rectMaxVol(matrix, options: RectMaxVolOptions(minRows: 3)) + + #expect(result.selectedRows == [0, 1, 3]) + #expect(result.iterations == 1) + #expect(result.converged) + try expectCoefficients( + result.coefficients, + rowMajorValues: [ + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.509_615_361_690_521_2, 0.221_153_840_422_630_3, 0.038_461_539_894_342_42, + 0.0, 0.0, 1.0, + ] + ) + try expectReconstruction(of: matrix, using: result) + try expectSelectedRowsAreIdentity(result) + } + + @Test func rectangularReferenceFixtureRespectsMaximumRowsLimit() throws { + let matrix = try DenseColumnMajorMatrix( + rows: 3, + columns: 2, + rowMajorValues: [ + 0.25, 0.0, + 0.5, 1.0, + 0.5, 1.5, + ] + ) + + let result = try rectMaxVol(matrix, options: RectMaxVolOptions(tolerance: 0.9, maxRows: 2)) + + #expect(result.selectedRows == [2, 0]) + #expect(result.iterations == 0) + #expect(!result.converged) + try expectCoefficients( + result.coefficients, + rowMajorValues: [ + 0.0, 1.0, + 0.666_666_626_930_236_8, 0.666_666_686_534_881_6, + 1.0, 0.0, + ] + ) + try expectReconstruction(of: matrix, using: result) + #expect(maximumUnselectedRowNorm(in: result) > 0.9) + } + + @Test func rectangularReferenceFixtureSelectsAllRowsWhenNeeded() throws { + let matrix = try DenseColumnMajorMatrix( + rows: 3, + columns: 2, + rowMajorValues: [ + 0.25, 0.0, + 0.5, 1.0, + 0.5, 1.5, + ] + ) + + let result = try rectMaxVol(matrix, options: RectMaxVolOptions(tolerance: 0.9, maxRows: 3)) + + #expect(result.selectedRows == [2, 0, 1]) + #expect(result.iterations == 1) + #expect(result.converged) + try expectCoefficients( + result.coefficients, + rowMajorValues: [ + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0, + 1.0, 0.0, 0.0, + ] + ) + try expectReconstruction(of: matrix, using: result) + try expectSelectedRowsAreIdentity(result) + } + + @Test func rectangularReferenceFixtureMatchesAppendCase() throws { + let matrix = try DenseColumnMajorMatrix( + rows: 5, + columns: 2, + rowMajorValues: [ + 1.0, 0.0, + 0.0, 1.0, + 1.2, 0.2, + 0.1, 1.3, + 0.8, 0.8, + ] + ) + + let result = try rectMaxVol(matrix, options: RectMaxVolOptions(minRows: 3)) + + #expect(result.selectedRows == [2, 3, 0]) + #expect(result.iterations == 1) + #expect(result.converged) + try expectCoefficients( + result.coefficients, + rowMajorValues: [ + 0.0, 0.0, 1.0, + 0.011_215_145_699_679_852, 0.767_505_407_333_374, -0.090_208_716_690_540_31, + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.399_453_848_600_387_6, 0.553_930_222_988_128_7, 0.265_262_335_538_864_14, + ] + ) + try expectReconstruction(of: matrix, using: result) + try expectSelectedRowsAreIdentity(result) + #expect(maximumUnselectedRowNorm(in: result) <= 1.0) + } + + @Test func randomizedOrthonormalTallMatrixSatisfiesFloatCriteria() throws { + let matrix = try orthonormalColumns(rows: 14, columns: 5, seed: 0xF10A7).mapValues(Float.init) + + let square = try maxVol(matrix, options: MaxVolOptions(tolerance: 1.05, maxIterations: 100)) + #expect(square.selectedRows.count == matrix.columns) + #expect(square.converged) + #expect(maximumAbsoluteCoefficient(in: square.coefficients) <= 1.05 + 1e-6) + try expectReconstruction(of: matrix, using: square, tolerance: 5e-5) + + let rectangular = try rectMaxVol(matrix, options: RectMaxVolOptions(tolerance: 1.0)) + #expect(rectangular.selectedRows.count >= matrix.columns) + #expect(rectangular.converged) + #expect(maximumUnselectedRowNorm(in: rectangular) <= 1.0 + 1e-5) + try expectReconstruction(of: matrix, using: rectangular, tolerance: 5e-5) + try expectSelectedRowsAreIdentity(rectangular) + } + + @Test func nearRankDeficientFloatInputThrowsDescriptiveError() throws { + let epsilon = Float.ulpOfOne / 16 + let matrix = try DenseColumnMajorMatrix( + rows: 3, + columns: 2, + rowMajorValues: [ + 1.0, 1.0, + 2.0, 2.0 + epsilon, + 3.0, 3.0, + ] + ) + + #expect(throws: MaxVolError.rankDeficient(pivot: 2)) { + try maxVol(matrix) + } + } +} + +private func expectReconstruction( + of matrix: DenseColumnMajorMatrix, + using result: MaxVolResult, + tolerance: Float = 1e-5 +) throws { + for row in 0.., + rowMajorValues: [Float], + tolerance: Float = 5e-6 +) throws { + let expected = try DenseColumnMajorMatrix( + rows: coefficients.rows, + columns: coefficients.columns, + rowMajorValues: rowMajorValues + ) + + for row in 0..) throws { + for (identityColumn, selectedRow) in result.selectedRows.enumerated() { + for column in 0.. +) -> Double { + coefficients.values.map { Double(abs($0)) }.max() ?? 0 +} + +private func maximumUnselectedRowNorm(in result: MaxVolResult) -> Double { + let selected = Set(result.selectedRows) + return (0.. DenseColumnMajorMatrix { + var generator = SeededGenerator(state: seed) + var columnVectors = (0.. [Double] in + (0.. Double { + zip(left, right).reduce(0) { total, pair in + total + pair.0 * pair.1 + } +} + +private struct SeededGenerator { + var state: UInt64 + + mutating func nextDouble() -> Double { + state = state &* 6_364_136_223_846_793_005 &+ 1_442_695_040_888_963_407 + let scaled = Double(state >> 11) / Double(UInt64.max >> 11) + return scaled * 2 - 1 + } +} + +private extension DenseColumnMajorMatrix where Scalar == Double { + func mapValues( + _ transform: (Double) -> Output + ) throws -> DenseColumnMajorMatrix { + try DenseColumnMajorMatrix( + rows: rows, + columns: columns, + columnMajorValues: values.map(transform) + ) + } +} diff --git a/Tests/MaxVolTests/NumericalBehaviorTests.swift b/Tests/MaxVolTests/NumericalBehaviorTests.swift new file mode 100644 index 0000000..bf9b587 --- /dev/null +++ b/Tests/MaxVolTests/NumericalBehaviorTests.swift @@ -0,0 +1,215 @@ +@testable import MaxVol +import Testing + +struct NumericalBehaviorTests { + @Test func randomizedOrthonormalMatricesSatisfyDoubleCriteriaAcrossSizes() throws { + for (rows, columns, seed) in [ + (8, 3, UInt64(0xD011_0001)), + (13, 4, UInt64(0xD011_0002)), + (21, 6, UInt64(0xD011_0003)), + ] { + let matrix = try orthonormalColumns(rows: rows, columns: columns, seed: seed) + + let square = try maxVol(matrix, options: MaxVolOptions(tolerance: 1.05, maxIterations: 200)) + #expect(square.selectedRows.count == matrix.columns) + #expect(square.converged) + #expect(maximumAbsoluteCoefficient(in: square.coefficients) <= 1.05) + try expectReconstruction(of: matrix, using: square, tolerance: 1e-9) + + let rectangular = try rectMaxVol(matrix, options: RectMaxVolOptions(tolerance: 1.0)) + #expect(rectangular.selectedRows.count >= matrix.columns) + #expect(rectangular.converged) + #expect(maximumUnselectedRowNorm(in: rectangular) <= 1.0 + 1e-12) + try expectReconstruction(of: matrix, using: rectangular, tolerance: 1e-9) + } + } + + @Test func randomizedOrthonormalMatricesSatisfyFloatCriteriaAcrossSizes() throws { + for (rows, columns, seed) in [ + (8, 3, UInt64(0xF10A_0001)), + (13, 4, UInt64(0xF10A_0002)), + (21, 6, UInt64(0xF10A_0003)), + ] { + let matrix = try orthonormalColumns(rows: rows, columns: columns, seed: seed).mapValues(Float.init) + + let square = try maxVol(matrix, options: MaxVolOptions(tolerance: 1.05, maxIterations: 200)) + #expect(square.selectedRows.count == matrix.columns) + #expect(square.converged) + #expect(maximumAbsoluteCoefficient(in: square.coefficients) <= 1.05 + 1e-6) + try expectReconstruction(of: matrix, using: square, tolerance: 5e-5) + + let rectangular = try rectMaxVol(matrix, options: RectMaxVolOptions(tolerance: 1.0)) + #expect(rectangular.selectedRows.count >= matrix.columns) + #expect(rectangular.converged) + #expect(maximumUnselectedRowNorm(in: rectangular) <= 1.0 + 1e-5) + try expectReconstruction(of: matrix, using: rectangular, tolerance: 5e-5) + } + } + + @Test func nearRankDeficientDoubleAndFloatInputsThrowDescriptiveErrors() throws { + let doubleMatrix = try DenseColumnMajorMatrix( + rows: 3, + columns: 2, + rowMajorValues: [ + 1.0, 1.0, + 2.0, 2.0 + Double.ulpOfOne / 16, + 3.0, 3.0, + ] + ) + let floatMatrix = try DenseColumnMajorMatrix( + rows: 3, + columns: 2, + rowMajorValues: [ + 1.0, 1.0, + 2.0, 2.0 + Float.ulpOfOne / 16, + 3.0, 3.0, + ] + ) + + #expect(throws: MaxVolError.rankDeficient(pivot: 2)) { + try maxVol(doubleMatrix) + } + #expect(throws: MaxVolError.rankDeficient(pivot: 2)) { + try maxVol(floatMatrix) + } + } +} + +private func expectReconstruction( + of matrix: DenseColumnMajorMatrix, + using result: MaxVolResult, + tolerance: Double +) throws { + for row in 0.., + using result: MaxVolResult, + tolerance: Float +) throws { + for row in 0.. +) -> Double { + coefficients.values.map(abs).max() ?? 0 +} + +private func maximumAbsoluteCoefficient( + in coefficients: DenseColumnMajorMatrix +) -> Double { + coefficients.values.map { Double(abs($0)) }.max() ?? 0 +} + +private func maximumUnselectedRowNorm(in result: MaxVolResult) -> Double { + let selected = Set(result.selectedRows) + return (0..) -> Double { + let selected = Set(result.selectedRows) + return (0.. DenseColumnMajorMatrix { + var generator = SeededGenerator(state: seed) + var columnVectors = (0.. [Double] in + (0.. Double { + zip(left, right).reduce(0) { total, pair in + total + pair.0 * pair.1 + } +} + +private struct SeededGenerator { + var state: UInt64 + + mutating func nextDouble() -> Double { + state = state &* 6_364_136_223_846_793_005 &+ 1_442_695_040_888_963_407 + let scaled = Double(state >> 11) / Double(UInt64.max >> 11) + return scaled * 2 - 1 + } +} + +private extension DenseColumnMajorMatrix where Scalar == Double { + func mapValues( + _ transform: (Double) -> Output + ) throws -> DenseColumnMajorMatrix { + try DenseColumnMajorMatrix( + rows: rows, + columns: columns, + columnMajorValues: values.map(transform) + ) + } +}