least_squares.jl 1.16 KB
Newer Older
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
39
40
41
42
43
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