Article

Solving Systems Using Iterative Methods

Use iterative methods to solve systems of equations where the coefficient matrix is sparse.

Overview

The code discussed here solves the following equation by using the iterative method Least Squares Minimum Residual (LSMR) to find the solution to this equation:

Equation to solve using the iterative method Least Squares Minimum Residual (LSMR).

In the equation above, A refers to the 4 x 3 matrix, and b to the right-hand-side vector. The code in this article solves the equation Ax = b by finding x.

Because A is overdetermined (it has more rows than columns), for most right-hand sides there will not be an exact solution. In this case, it is normally required to find the closest solution, in the sense that it minimizes the 2-norm of the error. That is, it solves the optimization min ‖ Ax - b ‖₂. This is known as the least-squares problem.

You could solve this problem through a direct method such as Sparse QR; however, for some problems, a faster method that provides a "good enough" solution is the iterative method LSMR. Unlike direct methods that factorize the matrix A, iterative methods only require the ability to multiply by the matrix (and its transpose, A). They move through a sequence of approximate solutions, converging to the correct answer. However, these methods run into numerical difficulties more often than direct methods. Resolving these issues requires expert knowledge, and is sometimes impossible. The most common method to improve and accelerate convergence is to use a preconditioner—an operator that approximates A¹, and is hence problem specific. For least-squares problems, just using a diagonal matrix with entries equal to the 2-norm of each column is often sufficient, and is the method discussed here.

Although you do not need the matrix explicitly to use iterative methods, the Accelerate framework offers the ability to supply a matrix explicitly if it is available, reducing the code required. The example takes advantage of this ability.

Create the Matrix

Use the code below—which is discussed in detail in Creating Sparse Matrices—to define the unsymmetric matrix A:

var rowIndices: [Int32] = [0, 1, 3,    // Column 0
                           0, 1, 2, 3, // Column 1
                           1, 2]       // Column 2
 
var values = [2.0, -0.2, 2.5,          // Column 0
              1.0, 3.2, -0.1, 1.1,     // Column 1
              1.4, 0.5]                // Column 2
 
var columnStarts = [0, 3, 7, 9]
 
let structure = SparseMatrixStructure(rowCount: 4,
                                      columnCount: 3,
                                      columnStarts: &columnStarts,
                                      rowIndices: &rowIndices,
                                      attributes: SparseAttributes_t(),
                                      blockSize: 1) 
let A = SparseMatrix_Double(structure: structure,
                            attributes: attributes,
                            blockSize: 1,
                            data: &values)

Solve the Equation

You define the right-hand-side and solution vectors as arrays wrapped in flexible data types. The initial values of x are used as an initial guess of the solution. If no good estimate is known, initialize all the values to zero:

var bValues = [1.200, 1.013, 0.205, -0.172]
let b = DenseVector_Double(count: 4,
                           data: &bValues) 
 
var xValues = [0.0, 0.0, 0.0]
let x = DenseVector_Double(count: 3,
                           data: &xValues)

With the matrix and vectors created, you are ready to perform the full LSMR iteration and iterate over the results in x. SparseSolve(_:_:_:_:_:) returns a status which, if equal to SparseIterativeConverged, indicates that the vector, x, contains the solution.

let status = SparseSolve(SparseLSMR(), A, b, x, SparsePreconditionerDiagScaling);
 
if status != SparseIterativeConverged {
    Swift.print("Failed to converge. Returned with error \(status)")
}
else {
    stride(from: 0, to: x.count, by: 1).forEach {i in
        Swift.print(x.data[Int(i)])
}}

You should get this output: x = 0.10 0.20 0.30.

See Also

Sparse Matrices

Creating Sparse Matrices

Create sparse matrices for factorization and solving systems.

Implementing Iterative Methods

Use iterative methods to solve large problems faster and with a lower memory overhead than with direct methods.

Solving Systems Using Direct Methods

Use direct methods to solve systems of equations where the coefficient matrix is sparse.

Creating a Sparse Matrix from Coordinate Format Arrays

Use separate coordinate format arrays to create sparse matrices.

Sparse Solvers

Solve systems of equations where the coefficient matrix is sparse.