From 6db679f28e22199dc25207270a66695663c32be4 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlsson <cont-frb@ulund.org> Date: Thu, 28 Jan 2016 15:48:17 +0100 Subject: [PATCH] implemented Tkhonov regularized least-squares using SVD --- src/least_squares.jl | 44 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 src/least_squares.jl diff --git a/src/least_squares.jl b/src/least_squares.jl new file mode 100644 index 0000000..23cf705 --- /dev/null +++ b/src/least_squares.jl @@ -0,0 +1,44 @@ +""" +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 + +function Lcurve(normE, normX, λ) + plot(normE,normX,xscale=:log10,yscale=:log10,m=:o) + annotations = [(normE[i],normX[i],"λ=$(round(λ[i],4))") for i in 1:length(λ)] + annotate!(annotations) + xlabel!("RMSE"); ylabel!("||k||"); title!("L-curve") +end -- GitLab