GitHub

Configuring Linear System Solvers

In the previous guide, we showed how to solve a consistent linear system in just a few lines of code. That example used the default configurations of the Kaczmarz solver solver, which is highly effective for many standard problems.

However, the true power of RLinearAlgebra.jl lies in its high degree of modularity and flexibility. You can fine-tune the solver's behavior by combining different ingredients, like cooking a fine dish, to tackle specific challenges, improve performance, or implement more complex algorithms.

We will follow the design philosophy of RLinearAlgebra.jl by composing different modules (Compressor, Logger, Solver, etc.) to customize a solver and solve the same consistent linear system,

\[Ax = b.\]

Configure Compressor

For large-scale problems, the matrix $A$ can be massive, slowing down iterative algorithms. We can use a randomized sketching technique to "compress" $A$ and $b$ to a lower dimension while preserving the essential information of the system, so that we can solve the system quickly (perhaps with some loss of accuracy).

Here, we'll configure a SparseSign compressor as an example. This compressor generates a sparse matrix $S$, whose non-zero elements are +1 or -1 (with scaling). The sparse matrix reduces the 1000 rows of our original system down to a more manageable $30$ rows.

# The goal is to compress the 1000 rows of A to 30 rows
compression_dim = 30
# We want each row of the compression matrix S to have 5 non-zero elements
non_zeros = 5

# Configure the SparseSign compressor
sparse_compressor = SparseSign(
    cardinality=Left(),                 # S will be left-multiplied: SAx = Sb
    compression_dim=compression_dim,    # The compressed dimension (number of rows)
    nnz=non_zeros,                      # The number of non-zero elements per row in S
    type=Float64                        # The element type for the compression matrix
)
SparseSign(Left(), 30, 5, Float64)

If the compression dimension of 30 rows is considered too large, it can be changed to 10 by updating the compressor configuration as follows:

# Change the dimension of the compressor. Similarly, you can use the same idea 
# for other configurations' changes.
sparse_compressor.compression_dim = 10;

The sparse_compressor contains all of the SparseSign ingredients needed to build a SparseSignRecipe.

(Optional) Build SparseSignRecipe and apply it

While the solver uses the sparse_compressor to construct the corresponding recipe, it is useful to construct the recipe at least once to see how it is used.

After defining the compressor's parameters, we combine it with our matrix $A$ to create a SparseSignRecipe. This "recipe" pre-calculates the sparse matrix S and prepares everything needed for efficient compression.

# Pass the compressor configuration and the original matrix A to
# create the final compression recipe.
S = complete_compressor(sparse_compressor, A);
SparseSignRecipe{Left}(Left(), 10, 1000, 5, [-0.4472135954999579, 0.4472135954999579], sparse([1, 4, 6, 7, 10, 2, 3, 4, 8, 9  …  2, 3, 4, 5, 8, 2, 4, 5, 6, 7], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2  …  999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000], [0.4472135954999579, 0.4472135954999579, 0.4472135954999579, 0.4472135954999579, -0.4472135954999579, 0.4472135954999579, 0.4472135954999579, 0.4472135954999579, -0.4472135954999579, -0.4472135954999579  …  -0.4472135954999579, -0.4472135954999579, -0.4472135954999579, 0.4472135954999579, -0.4472135954999579, -0.4472135954999579, -0.4472135954999579, -0.4472135954999579, 0.4472135954999579, -0.4472135954999579], 10, 1000))

We can then apply S to matrices and vectors through multiplication. Other multiplication functionality is available as well (i.e., mul!).

# Form the compressed system SAx = Sb
SA = S * A;
Sb = S * b;
10-element Vector{Float64}:
   3.3440916591932424
  -3.063659176453057
  -7.4403916470616185
 -25.590506108373184
 -57.77254328823237
  23.79198530341466
  37.0239963105425
  56.09803546431352
  -4.917934174810977
 -14.725095900001218

Configure Logger

To monitor the solver and control its execution, we will configure a BasicLogger . This object serves two purposes: tracking metrics (like the error history) and defining stopping rules.

We'll configure it to stop after a maximum of 500 iterations or if the residual error drops below 1e-6. We will also set collection_rate = 5 to record the error every $5$ iterations.

# Configure the logger to control the solver's execution
logger = BasicLogger(
    max_it = 500,
    threshold = 1e-6,
    collection_rate = 5
)
BasicLogger(500, 5, 1.0e-6, RLinearAlgebra.threshold_stop)

(Optional) Build the BasicLoggerRecipe

Just like the compressor, this logger is just a set of instructions. To make it "ready" to store data, we could call complete_logger on it:

# We can create the recipe manually, though this is rarely needed
logger_recipe = complete_logger(logger)
BasicLoggerRecipe{typeof(threshold_stop)}(500, 0.0, 1.0e-6, 1, 1, 5, false, RLinearAlgebra.threshold_stop, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

The logger_recipe's has fields for storing the error history, the current iteration, and the converged status.

Again, you almost never need to call complete_logger yourself. Because the solver (which we will configure next) or the rsolve! (the function that solves the system) handles it for us. When we call complete_solver or the rsolve!, it will automatically find the logger inside the solver, call complete_logger on it, and store the resulting logger_recipe inside the final solver_recipe.

Configure Solver

Now we assemble our configured ingredients—the sparse_compressor and the logger—into the main Kaczmarz solver object. For any component we don't specify, a default will be used. Here, we'll explicitly specify the LQSolver as our sub-solver.

# Create the Kaczmarz solver object by passing in the ingredients
kaczmarz_solver = Kaczmarz(
    compressor = sparse_compressor,
    log = logger,
    sub_solver = LQSolver()
)
Kaczmarz(1.0, SparseSign(Left(), 10, 5, Float64), BasicLogger(500, 5, 1.0e-6, RLinearAlgebra.threshold_stop), FullResidual(), LQSolver())

(Optional) create the SolverRecipe

This is the step where everything comes together. Just as with the compressor, when we call complete_solver, it takes the kaczmarz_solver configurations and the problem data, and then:

  1. Pre-allocates memory for everything needed in the algorithm.

  2. Finds the sparse_compressor config inside kaczmarz_solver and calls complete_compressor to create the SparseSignRecipe.

  3. Finds the logger inside kaczmarz_solver and calls complete_logger to create the BasicLoggerRecipe.

  4. Bundles all these "recipes" into a single, ready-to-use solver_recipe.

# Set the solution vector x (typically a zero vector)
solution = zeros(Float64, num_cols);
# Create the solver recipe by combining the solver and the problem data
solver_recipe = complete_solver(kaczmarz_solver, solution, A, b);
KaczmarzRecipe{Float64, Vector{Float64}, Matrix{Float64}, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, false}, SparseSignRecipe{Left}, BasicLoggerRecipe{typeof(threshold_stop)}, FullResidualRecipe{Vector{Float64}}, LQSolverRecipe{Matrix{Float64}}}(SparseSignRecipe{Left}(Left(), 10, 1000, 5, [-0.4472135954999579, 0.4472135954999579], sparse([4, 6, 8, 9, 10, 1, 4, 6, 7, 10  …  2, 4, 5, 9, 10, 2, 3, 4, 6, 10], [1, 1, 1, 1, 1, 2, 2, 2, 2, 2  …  999, 999, 999, 999, 999, 1000, 1000, 1000, 1000, 1000], [0.4472135954999579, -0.4472135954999579, 0.4472135954999579, 0.4472135954999579, 0.4472135954999579, 0.4472135954999579, -0.4472135954999579, -0.4472135954999579, -0.4472135954999579, 0.4472135954999579  …  -0.4472135954999579, 0.4472135954999579, -0.4472135954999579, 0.4472135954999579, 0.4472135954999579, 0.4472135954999579, -0.4472135954999579, 0.4472135954999579, -0.4472135954999579, 0.4472135954999579], 10, 1000)), BasicLoggerRecipe{typeof(threshold_stop)}(500, 0.0, 1.0e-6, 1, 1, 5, false, RLinearAlgebra.threshold_stop, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), FullResidualRecipe{Vector{Float64}}([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), LQSolverRecipe{Matrix{Float64}}([0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]), 1.0, [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

Again, this step can also be skipped and directly pass the solver config kaczmarz_solver into the function that can solve the system.

Solve and Verify the Result

Solve the System

We call rsolve! to run the Kaczmarz solver. The ! in the name indicates that the function modifies its arguments in-place. Here, solution will be updated with the solution vector. The algorithm will run until a stopping criterion from our logger is met.

# Set the solution vector x (typically a zero vector)
solution = zeros(Float64, num_cols);

# Run the solver!
_, solver_history = rsolve!(kaczmarz_solver, solution, A, b);

# The solution is now stored in the updated solution vector
solution
50-element Vector{Float64}:
  0.7616097937555282
 -1.0226412740402329
 -1.7481863019033341
 -1.8614360313665081
  0.7521297781149382
 -0.8366637730361783
 -1.7506683920011195
  0.9523824948544759
 -0.4725412995357882
 -0.869265815497822
  ⋮
 -1.0914581593342008
 -1.1900315749895503
  0.27053311915010925
 -1.4376970577578472
 -0.25589560651948495
 -0.3674764808565056
  0.5623435372178435
 -0.49560261223091634
 -0.8351135428945264

Verify the result

Finally, let's check how close our calculated solution is to the known x_true and inspect the logger to see how the solver performed.

# We can inspect the logger's history to see the convergence
error_history = solver_history.log.hist;
println(" - Solver stopped at iteration: ", solver_history.log.iteration)
println(" - Final residual error, ||Ax-b||_2: ", error_history[solver_history.log.record_location])

# Calculate the norm of the error
error_norm = norm(solution - x_true)
println(" - Norm of the error between the solution and x_true: ", error_norm)
 - Solver stopped at iteration: 175
 - Final residual error, ||Ax-b||_2: 8.641657903476489e-7
 - Norm of the error between the solution and x_true: 2.7693514309706242e-8

As you can see, by composing different kernels, we configured a randomized solver that found a solution vector very close to the true solution.