A new package: RandomLinearAlgebraSolvers.jl

1 minute read

Published:

The package RandomLinearAlgebraSolvers.jl is now a Julia package. It contains randomized iterative methods for linear algebra and random projectors used to solve linear systems in the sense of the Johnson–Lindenstrauss lemma. The package is still at an early stage and new contributions are very welcome. You will find below an example extracted from the package’s documentation. The package uses Stopping.jl as a framework for iterative methods.

Benchmark on dense overdetermined random matrices

Benchmark of methods to solve Ax = b with A a randomly generated m x n matrix, and b = A * xref.

using RandomLinearAlgebraSolvers, Random, LinearAlgebra, SparseArrays
Random.seed!(1234)
using DataFrames, Printf, SolverBenchmark, Stopping

N = 5 # number of problems
m, n = 10000, 100 # size of A: m x n

names = [:RandomizedKaczmarz, :RandomizedBlockKaczmarz, :RandomizedCD]

#Initialization of the DataFrame for n problems.
stats = Dict(name => DataFrame(
         :id     => 1:N,
         :nvar   => zeros(Int64, N),
         :status => [:Unknown for i = 1:N],
         :time   => NaN*ones(N),
         :iter   => zeros(Int64, N),
         :score  => NaN*ones(N)) for name in names)

for i=1:N

  A = 100 * rand(m, n)
  xref = 100 * rand(n)
  b = A * xref

  x0 = zeros(size(A,2))
  la_stop = LAStopping(A, b, max_iter = 100000, rtol = sqrt(eps()), atol = sqrt(eps()))
  for name in names

    #solve the problem
    reinit!(la_stop, rstate = true, x = x0, res = similar(b))
    la_stop.meta.start_time = time()
    @time eval(name)(la_stop, r = 80, is_zero_start=true)
    sol = la_stop.current_state.x

    #update the stats from the Stopping
    stats[name].nvar[i]   = n
    stats[name].status[i] = status(la_stop)
    stop_has_time = (la_stop.current_state.current_time != nothing)
    stats[name].time[i]   =  stop_has_time ? la_stop.current_state.current_time - la_stop.meta.start_time : time() - la_stop.meta.start_time
    stats[name].iter[i]   = la_stop.meta.nb_of_stop
    stats[name].score[i]  = norm(la_stop.current_state.current_score, Inf)

  end

end

using Plots
gr()
cost(df) = (df.status .!= :Optimal) * Inf + df.time
p = performance_profile(stats, cost)