Dorr Matrix

metnum_url = "https://github.com/santiago-imelio/metnum.git"

Mix.install([
  {:metnum, git: metnum_url, branch: "main"},
  {:exla, "~> 0.7.3"},
  {:tucan, "~> 0.3.1"},
  {:kino_vega_lite, "~> 0.1.13"}
])

Nx.global_default_backend(EXLA.Backend)
Nx.Defn.global_default_options(compiler: EXLA)

Introduction

In this notebook, we test Jacobi, Gauss-Seidel and SOR iterative methods against the Dorr matrix, a diagonally dominant, tridiagonal, M-matrix. Then, we will compare the three methods by visualizing the evolution of the residual and step error.

dorr = Metnum.Matrix.dorr(5, 1.0)

We define a module with some utility functions to plot the evolution of the iterative methods.

alias Metnum.LinearEquations

defmodule Plot do
  def error_evolution(seq, solver, epochs, opts \\ []) do
    errors =
      Nx.subtract(seq[1..(epochs - 1)], seq[0..(epochs - 2)])
      |> Nx.vectorize(:deltas)
      |> Nx.LinAlg.norm(ord: 2)
      |> Nx.devectorize(keep_names: false)
      |> Nx.to_list()

    default_opts = [
      title: "Evolution of iteration error",
      height: 400,
      width: 700,
      line_color: "green"
    ]

    opts = Keyword.merge(default_opts, opts)

    data_errors = [
      error: errors,
      epoch: 0..(epochs - 1),
      solver: List.duplicate(solver, epochs)
    ]

    Tucan.lineplot(data_errors, "epoch", "error", opts)
  end

  def residual_evolution(seq, solver, solution, epochs, opts \\ []) do
    residuals =
      Nx.subtract(seq[1..(epochs - 1)], solution)
      |> Nx.vectorize(:deltas)
      |> Nx.LinAlg.norm(ord: 2)
      |> Nx.devectorize(keep_names: false)
      |> Nx.to_list()

    default_opts = [
      title: "Evolution of residual",
      height: 400,
      width: 700,
      line_color: "orange"
    ]

    opts = Keyword.merge(default_opts, opts)

    data_residuals = [
      residual: residuals,
      epoch: 0..(epochs - 1),
      solver: List.duplicate(solver, epochs)
    ]

    Tucan.lineplot(data_residuals, "epoch", "residual", opts)
  end
end

We define a Dorr matrix of size 500 by 500 and $\theta = 3$, which we will use to run our iterative methods against.

a = Metnum.Matrix.dorr(500, 3)

Here we define the parameters that we will use for Jacobi, Gauss-Seidel and SOR. Experiments will hav a maximum of 50,000 iterations per run and a tolerance of 0.001 for the difference between the previous and current value.

b = Nx.iota({500})
x0 = Nx.tile(Nx.tensor(0.0), [500])
opts = [max_epochs: 50_000, tolerance: 0.001, sequence: true]

We will use the real solution to the problem $Ax = b$ to see the evolution of the residual error.

solution = Nx.LinAlg.solve(a, b)

Running experiments concurrently

We run each iterative method using the previous parameters. Since each experiment will take on its own about 40 seconds, we’ll speed things up by taking advantage of Elixir’s concurrency, using the Task module.

To do this we build a list of all the solvers and pass it to Task.async_stream/3. This function will spawn a process per enum item to evaluate the given fun concurrently. In this case, the passed fun will update opts accordingly and run the iterative method with the given solver.

Finally we apply Enum.map/2 to execute the async stream and map the results.

solvers = [:jacobi, :gauss_seidel, :sor]

results = 
  solvers
  |> Task.async_stream(fn solver ->
    opts = Keyword.put(opts, :solver, solver)
    {solver, LinearEquations.solve(a, b, x0, opts)}
  end, 
  [timeout: :infinity]
)
|> Enum.map(fn {:ok, result} -> result end)

As we see, for each iterative method we returned a tuple with the solver and the experiment results. This is because the order of the results from async_stream/3 is nondeterministic, and we want to know which results correspond to each solver.

{_, {jacobi_seq, _}} = Enum.find(results, fn {solver, _} -> solver == :jacobi end)
{_, {gs_seq, _}} = Enum.find(results, fn {solver, _} -> solver == :gauss_seidel end)
{_, {sor_seq, _}} = Enum.find(results, fn {solver, _} -> solver == :sor end)

Comparing solvers

Tucan.vconcat([
  Tucan.layers([
    Plot.error_evolution(jacobi_seq, "Jacobi", opts[:max_epochs]),
    Plot.error_evolution(gs_seq, "Gauss-Seidel", opts[:max_epochs]),
    Plot.error_evolution(sor_seq, "SOR", opts[:max_epochs])
  ])
  |> Tucan.color_by("solver"),
  Tucan.layers([
    Plot.residual_evolution(jacobi_seq, "Jacobi", solution, opts[:max_epochs]),
    Plot.residual_evolution(gs_seq, "Gauss-Seidel", solution, opts[:max_epochs]),
    Plot.residual_evolution(sor_seq, "SOR", solution, opts[:max_epochs])
  ])
  |> Tucan.color_by("solver")
])
|> Tucan.set_theme(:google_charts)