least_squares.jl 1.25 KB
 Fredrik Bagge Carlson committed Jan 28, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 ``````""" 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 `````` Fredrik Bagge Carlson committed Sep 13, 2016 39 40 41 42 43 ``````""" Lcurve(normE, normX, λ) Plots the L-curve. normE and normX are obtained from, e.g., `dsvd` """ `````` Fredrik Bagge Carlson committed Jan 28, 2016 44 ``````function Lcurve(normE, normX, λ) `````` Fredrik Bagge Carlson committed Sep 13, 2016 45 `````` plot(normE,normX,xscale=:log10,yscale=:log10,m=:o, xlabel="RMSE", ylabel="||k||", title="L-curve") `````` Fredrik Bagge Carlson committed Sep 13, 2016 46 `````` annotations = [(normE[i],normX[i],"λ=\$(round(λ[i],8))") for i in 1:length(λ)] `````` Fredrik Bagge Carlson committed Jan 28, 2016 47 48 `````` annotate!(annotations) end``````