Found the solution... The problem was the Hashable conformance. After removing it and changing ContentView so that it can work without Hashable models, the problem was gone.

Found an easier solution: using LAPACK's dgels function. Thanks for your help anyway, @Claude31 !

Running that code works fine with the matrix they used in the example, but as soon as I try it out with other matrices, the function crashes with this error: 'Matrix is structurally singular'. Which, according to Google, means that the matrix is not invertible. So seemingly QR-factorization (the factorization used in the example) only works with matrices that are invertible. I tried out the other methods of factorization, but they all have their limitations; none of them seem suitable for what I'm trying to achieve. However, the other leastSquaresSolution-method that I initially wrote also returns zero's if I initialize it with nan's instead of zero's. So xPtr must somehow be modified by the SparseSolve function (although I also have no idea how that's possible). Moreover, exactly the same code (literally copy-pasted) returns the right values when I call it in my real project, except that it sometimes throws an EXC_BAD_ACCESS. This stuff is far more complicated than I had expected it to be :)

Thank you for your reaction! The code you mentioned is derived from this article, where the xPtr isn’t changed manually. If I understood it correctly, the SparseSolve function changes x, which in turn changes xPtr and xValues.

the sparseMatrix function with comments: extension Array where Element == [Double] { // A sparse matrix is a mattrix where all zero's are ommitted. // Normal matrix: Sparse matrix: // 0 1 0 1 // 1 3 0 1 3 // 4 0 2 4 2 // Find the sparse matrix of this matrix. Returns an array of doubles containing the values of the sparse matrix (omitting all zero's) and a SparseMatrixStructure containing information about where the columns start and where the rows start. The values array of the sparse matrix above would be [1, 4, 1, 3, 2]. func sparseMatrix() -> (structure: SparseMatrixStructure, values: [Double]) { let columns = self.transpose() // Get the row indices of the matrix. The row indices of the sparse matrix above would be // [1, 2 column 0 // 0, 1, column 1 // 2] column 2 var rowIndices: [Int32] = columns.map { column in column.indices.compactMap { indexInColumn in if column[indexInColumn] != 0 { return Int32(indexInColumn) } return nil } }.reduce([], +) // Get the column starts. For the sparse matrix above, the column starts would be // [1, column 0 // 3, column 1 // 5, column 2 // 5] the total amount of values in the sparse matrix let sparseColumns = columns.map { column in column.compactMap { if $0 != 0 { return $0 } return nil } } var counter = 0 var columnStarts = [Int]() for sparseColumn in sparseColumns { columnStarts.append(counter) counter += sparseColumn.count } // Get the reducedSparseColumns (the `values`). For the matrix above, these would be [1, 4, 1, 3, 2]. let reducedSparseColumns = sparseColumns.reduce([], +) columnStarts.append(reducedSparseColumns.count) // Wrap the information about the sparse matrix in a SparseMatrixStructure let structure: SparseMatrixStructure = rowIndices.withUnsafeMutableBufferPointer { rowIndicesPtr in columnStarts.withUnsafeMutableBufferPointer { columnStartsPtr in let attributes = SparseAttributes_t() return SparseMatrixStructure( rowCount: Int32(self.count), columnCount: Int32(columns.count), columnStarts: columnStartsPtr.baseAddress!, rowIndices: rowIndicesPtr.baseAddress!, attributes: attributes, blockSize: 1 ) } } // Return the structure and the values of the sparse matrix. return (structure, reducedSparseColumns) } // swaps rows and columns func transpose() -> Self { let columns = self.count let rows = self.reduce(0) { Swift.max($0, $1.count) } return (0 ..< rows).reduce(into: []) { result, row in result.append((0 ..< columns).reduce(into: []) { result, column in result.append(row < self[column].count ? self[column][row] : 0) }) } } }

The SavitzkyGolay class with comments: class SavitzkyGolay { static func coefficients(windowLength: Int, polynomialOrder: Int, derivativeOrder: Int = 0, delta: Int = 1) -> [Double] { // windowLength is the number of coefficients returned by this function. guard windowLength > 0 else { fatalError("windowLength must be positive") } // polynomialOrder is the order of the polynomial used to smooth the noisy data (the coefficients are calculated independently from the noisy data) guard polynomialOrder < windowLength else { fatalError("polynomialOrder must be less than windowLength") } // The derivativeOrder is the order of the derivative to compute. If it is set to zero, the noisy data is smoothed without differentiating. guard derivativeOrder <= polynomialOrder else { return [Double](repeating: 0, count: windowLength) } let (halfWindow, remainder) = windowLength.quotientAndRemainder(dividingBy: 2) var pos = Double(halfWindow) // pos should not be a round number (because otherwise this function won't work appropriately), so if windowLength is even, which means that halfWindow and thereby pos are round numbers, subtract 0.5 from pos. if remainder == 0 { pos -= 0.5 } // Form the matrix A. The columns of A are powers of the integers from -pos to windowLength - pos - 1. The powers (i.e., rows) range from 0 to polynomialOrder. let V = [Double](stride(from: Double(windowLength) - pos - 1, through: -pos, by: -1)) let P = [Double](stride(from: 0, through: Double(polynomialOrder), by: 1)) let A = P.map { exponent in V.map { pow($0, exponent) } } // Form the vector B. It determines which order derivative is returned. var B = [Double](repeating: 0, count: polynomialOrder + 1) // The value assigned to B[derivativeOrder] scales the result to take into account the order of the derivative and the distance between two values of the noisy data. B[derivativeOrder] = Double(factorial(derivativeOrder)) / pow(Double(delta), Double(derivativeOrder)) // Return the least-squares solution of Ax = B return leastSquaresSolution(A: A, B: B) } } leastSquaresSolution with comments: func leastSquaresSolution(A: [[Double]], B: [Double]) -> [Double] { // Get the sparse matrix of A. A sparse matrix is a matrix that omits all zero's. sparseMatrix() returns an array of doubles, omitting all zero's, together with information about where in the array the columns start and what the indices of the rows are (for more information, take a look at the sparseMatrix() function). let sparseA = A.sparseMatrix() // Copy some constants and save them in variables, so that they can be pointed to using withUnsafeMutableBufferPointer. var sparseAValuesCopy = sparseA.values var bValues = B // The vector X returned by this function has as many values as A has columns. In other words, X has as many values as A's transpose has rows. So X should be initialised to have that amount of values. var xValues = [Double](repeating: 0, count: A.transpose().count) sparseAValuesCopy.withUnsafeMutableBufferPointer { valuesPtr in // Construct the matrix A in Ax = B let a = SparseMatrix_Double( structure: sparseA.structure, data: valuesPtr.baseAddress! ) bValues.withUnsafeMutableBufferPointer { bPtr in xValues.withUnsafeMutableBufferPointer { xPtr in // Construct the matrices B and x in Ax = B let b = DenseVector_Double( count: Int32(B.count), data: bPtr.baseAddress! ) let x = DenseVector_Double( count: Int32(A.transpose().count), data: xPtr.baseAddress! ) // Try to find X with SparseSolve (this is where it seems to go wrong) let status = SparseSolve(SparseLSMR(), a, b, x, SparsePreconditionerDiagScaling) if status != SparseIterativeConverged { fatalError("Failed to converge. Returned with error \(status).") } } } } // Should return x in Ax = B return xValues }

Hello Claude, First of all, I am sorry about the absence of comments. It would be a headache for me too if I had not written the code myself, so I should have thought about adding comments earlier... I will add them and send the code again soon. The xValues form the vector x in the equation Ax = B, where A and B are matrices that are the arguments of the leastSquaresSolution function. Apple describes how the Ax = B equation can be solved with Swift in this article. xValues is indeed modified by updating xPtr. xPtr's baseAddress is passed to the DenseVector_Double that holds x: let x = DenseVector_Double( count: Int32(A.transpose().count), data: xPtr.baseAddress! ) Then, x is passed to the constant 'status', which should modify x and thereby should modify xPtr and xValues: let status = SparseSolve(SparseLSMR(), a, b, x, SparsePreconditionerDiagScaling)

Things get even more interesting, as I just found out that the leastSquaresSolution function always returns an array of zero's, instead of the values it should return, while it does actually return the right values in my real code...