using Krylov
using LinearAlgebra, Printf
m = 5
n = 8
λ = 1.0e-3
A = rand(m, n)
b = A * ones(n)
xy_exact = [A λ*I] \ b # In Julia, this is the min-norm solution!
(x, y, stats) = craig(A, b, λ=λ, atol=0.0, rtol=1.0e-20, verbose=1)
show(stats)
# Check that we have a minimum-norm solution.
# When λ > 0 we solve min ‖(x,s)‖ s.t. Ax + λs = b, and we get s = λy.
@printf("Primal feasibility: %7.1e\n", norm(b - A * x - λ^2 * y) / norm(b))
@printf("Dual feasibility: %7.1e\n", norm(x - A' * y) / norm(x))
@printf("Error in x: %7.1e\n", norm(x - xy_exact[1:n]) / norm(xy_exact[1:n]))
if λ > 0.0
@printf("Error in y: %7.1e\n", norm(λ * y - xy_exact[n+1:n+m]) / norm(xy_exact[n+1:n+m]))
end
CRAIG: system of 5 equations in 8 variables
k ‖r‖ ‖x‖ ‖A‖ κ(A) α β timer
0 1.01e+01 0.00e+00 0.00e+00 0.00e+00 ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ 0.38s
1 3.55e-01 2.78e+00 3.63e+00 3.63e+00 3.6e+00 1.3e-01 0.79s
2 4.98e-02 2.80e+00 3.82e+00 5.41e+00 1.2e+00 1.7e-01 0.93s
3 5.78e-03 2.80e+00 3.90e+00 6.76e+00 7.5e-01 8.8e-02 0.93s
4 2.42e-03 2.80e+00 3.92e+00 7.86e+00 4.0e-01 1.7e-01 0.93s
5 8.65e-10 2.80e+00 3.99e+00 9.01e+00 7.2e-01 2.6e-07 0.93s
SimpleStats
niter: 5
solved: true
inconsistent: false
residuals: []
Aresiduals: []
κ₂(A): []
timer: 932.47ms
status: solution good enough for the tolerances given
Primal feasibility: 8.6e-11
Dual feasibility: 2.6e-16
Error in x: 8.5e-11
Error in y: 8.0e-11