IRLS

using Plots
using LinearAlgebra: norm, opnorm, pinv
using Interact
using Random: seed!

The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.

m, n = 100, 50
seed!(0)
A = randn(m,n)
b = randn(m)
σ1 = opnorm(A); # largest singular value is the operator 2-norm of A

xls = pinv(A)*b # SVD based solution
iters = 1000

@manipulate for μ_factor in range(0.0, 2.5, 51)

    μ = μ_factor / σ1^2

    x = zeros(size(A,2))
    ils_err = zeros(iters+1)

    ils_err[1] = norm(x - xls)
    for k = 1:iters
        x = x - μ * (A'*(A*x-b))
        ils_err[k+1] = norm(x - xls) / norm(xls)
    end
    scatter(0:iters, log10.(ils_err),
        xlabel="Iteration", xtick=[0,iters],
        ylabel="log10(NRMSD)", ylim=(-15,5),
        label = "μ_factor = $(round(μ_factor,digits=3))",
        title="Iterative LS convergence plot")

end

plot!() # for saving as html