A new package: RandomLinearAlgebraSolvers.jl
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)
