least_squares.jl

"""
Calculates the Tikhonov regularized estimate min. ||Ax-b|| + λ||x|| using SVD
for a vector with different values of λ

`x, normE, normX = dsvd(A,U,S,V,b,λ)`

`x, normE, normX = dsvd(A,b,λ)`
"""
function dsvd(A,U,S,V,b,λ)
    if minimum(λ) < 0
        error("Illegal regularization parameter λ")
    end
    UTb     = U'b
    N       = size(U,1)
    n       = size(V,1)
    nl      = length(λ)
    x       = zeros(n,nl)
    normE   = zeros(nl)
    normX   = zeros(nl)
    print("dsvd: ")
    for (i,λi) in enumerate(λ)
        x[:,i] = V*diagm(S./(S.^2+λi))*UTb
        normX[i] = norm(x[:,i])
        normE[i] = norm(A*x[:,i] - b)
        print("[$(round(i/length(λ),2)),$(round(λi,2))]  ")
    end
    println("Done ")
    return x, normE, normX
end

function dsvd(A,b,λ)
    if minimum(λ) < 0
        error("Illegal regularization parameter λ")
    end
    U,S,V = svd(A)
    dsvd(A,U,S,V,b,λ)
end

"""
    Lcurve(normE, normX, λ)

Plots the L-curve. normE and normX are obtained from, e.g., `dsvd`
"""
function Lcurve(normE, normX, λ)
    plot(normE,normX,xscale=:log10,yscale=:log10,m=:o, xlabel="RMSE", ylabel="||k||", title="L-curve")
    annotations = [(normE[i],normX[i],"λ=$(round(λ[i],8))") for i in 1:length(λ)]
    annotate!(annotations)
end