[Master Index] [Index for PublicToolbox/regutools]

# regudemo

## (PublicToolbox/regutools/regudemo.m in BrainStorm 2.0 (Alpha))

### Help Text

    REGUDEMO Tutorial script for Regularization Tools.

### Cross-Reference Information

This calls

    csvd                               C:\BrainStorm_2001\PublicToolbox\regutools\csvd.m
    fil_fac                            C:\BrainStorm_2001\PublicToolbox\regutools\fil_fac.m
    gcv                                C:\BrainStorm_2001\PublicToolbox\regutools\gcv.m
    get_l                              C:\BrainStorm_2001\PublicToolbox\regutools\get_l.m
    gsvd                               C:\BrainStorm_2001\PublicToolbox\regutools\gsvd.m
    ilaplace                           C:\BrainStorm_2001\PublicToolbox\regutools\ilaplace.m
    l_curve                            C:\BrainStorm_2001\PublicToolbox\regutools\l_curve.m
    lsqr                               C:\BrainStorm_2001\PublicToolbox\regutools\lsqr.m
    picard                             C:\BrainStorm_2001\PublicToolbox\regutools\picard.m
    plot_lc                            C:\BrainStorm_2001\PublicToolbox\regutools\plot_lc.m
    shaw                               C:\BrainStorm_2001\PublicToolbox\regutools\shaw.m
    tgsvd                              C:\BrainStorm_2001\PublicToolbox\regutools\tgsvd.m
    tikhonov                           C:\BrainStorm_2001\PublicToolbox\regutools\tikhonov.m
    tsvd                               C:\BrainStorm_2001\PublicToolbox\regutools\tsvd.m
    ursell                             C:\BrainStorm_2001\PublicToolbox\regutools\ursell.m

### Listing of C:\BrainStorm_2001\PublicToolbox\regutools\regudemo.m

    %REGUDEMO Tutorial script for Regularization Tools.

    % Per Christian Hansen, UNI-C, 03/17/93.

    echo on, clf

    % Part 1.  The discrete Picard condition
    % --------------------------------------
    %
    % First generate a "pure" test problem where only rounding
    % errors are present.  Then generate another "noisy" test
    % problem by adding white noise to the right-hand side.
    %
    % Next compute the SVD of the coefficient matrix A.
    %
    % Finally, check the Picard condition for both test problems
    % graphically.  Notice that for both problems the condition is
    % indeed satisfied for the coefficients corresponding to the
    % larger singular values, while the noise eventually starts to
    % dominate.

    [A,b_bar,x] = shaw(32);
    randn('seed',70957);
    e = 1e-3*randn(size(b_bar)); b = b_bar + e;
    [U,s,V] = csvd(A);
    pause % Strike any key to continue
    subplot(2,1,1); picard(U,s,b_bar);
    subplot(2,1,2); picard(U,s,b); pause
    clf

    % Part 2.  Filter factors
    % -----------------------
    %
    % Compute regularized solutions to the "noisy" problem from Part 1
    % by means of Tikhonov's method and LSQR without reorthogonalization.
    % Also, compute the corresponding filter factors.
    %
    % A surface (or mesh) plot of the solutions clearly shows their dependence
    % on the regularization parameter (lambda or the iteration number).

    lambda = [1,3e-1,1e-1,3e-2,1e-2,3e-3,1e-3,3e-4,1e-4,3e-5];
    X_tikh = tikhonov(U,s,V,b,lambda);
    F_tikh = fil_fac(s,lambda);
    iter = 30; reorth = 0;
    [X_lsqr,rho,eta,F_lsqr] = lsqr(A,b,iter,reorth,s);
    pause % Strike any key to continue
    subplot(2,2,1); surf(X_tikh), title('Tikhonov solutions'), axis('ij')
    subplot(2,2,2); surf(log10(F_tikh)), axis('ij')
                    title('Tikh filter factors, log scale')
    subplot(2,2,3); surf(X_lsqr(:,1:17)), title('LSQR solutions'), axis('ij')
    subplot(2,2,4); surf(log10(F_lsqr(:,1:17))), axis('ij')
                    title('LSQR filter factors, log scale'), pause
    clf

    % Part 3.  The L-curve
    % --------------------
    %
    % Plot the L-curves for Tikhonov regularization and for
    % LSQR for the "noisy" test problem from Part 1.
    %
    % Notice the similarity between the two L-curves and thus,
    % in turn, by the two methods.

    pause % Strike any key to continue
    subplot(1,2,1); l_curve(U,s,b); axis([1e-3,1,1,1e3])
    subplot(1,2,2); plot_lc(rho,eta,'o'); axis([1e-3,1,1,1e3]), pause
    clf

    % Part 4.  Regularization parameters
    % ----------------------------------
    %
    % Use the L-curve criterion and GCV to determine the regularization
    % parameters for Tikhonov regularization and truncated SVD.
    %
    % Then compute the relative errors for the four solutions.

    pause % Strike any key to continue
    lambda_l = l_curve(U,s,b);   axis([1e-3,1,1,1e3]),      pause
    k_l = l_curve(U,s,b,'tsvd'); axis([1e-3,1,1,1e3]),      pause
    lambda_gcv = gcv(U,s,b);     axis([1e-6,1,1e-9,1e-1]),  pause
    k_gcv = gcv(U,s,b,'tsvd');   axis([0,20,1e-9,1e-1]),    pause

    x_tikh_l   = tikhonov(U,s,V,b,lambda_l);
    x_tikh_gcv = tikhonov(U,s,V,b,lambda_gcv);
    if isnan(k_l)
      x_tsvd_l = zeros(32,1); % Spline Toolbox not available.
    else
      x_tsvd_l = tsvd(U,s,V,b,k_l);
    end
    x_tsvd_gcv = tsvd(U,s,V,b,k_gcv);
    [norm(x-x_tikh_l),norm(x-x_tikh_gcv),...
     norm(x-x_tsvd_l),norm(x-x_tsvd_gcv)]/norm(x)
    pause % Strike any key to continue

    % Part 5.  Standard form versus general form
    % ------------------------------------------
    %
    % Generate a new test problem: inverse Laplace transformation
    % with white noise in the right-hand side.
    %
    % For the general-form regularization, choose minimization of
    % the first derivative.
    %
    % First display some left singular vectors of SVD and GSVD; then
    % compare truncated SVD solutions with truncated GSVD solutions.
    % Notice that TSVD cannot reproduce the asymptotic part of the
    % solution in the right part of the figure.

    n = 16; [A,b,x] = ilaplace(n,2);
    b = b + 1e-4*randn(size(b));
    L = get_l(n,1);
    [U,s,V] = csvd(A); [UU,VV,sm,XX] = gsvd(A,L);
    pause % Strike any key to continue
    I = 1;
    for i=[3,6,9,12]
      subplot(2,2,I); plot(1:n,V(:,i)); axis([1,n,-1,1])
      xlabel(['i = ',num2str(i)]), I = I + 1;
    end
    subplot(2,2,1), text(12,1.2,'Right singular vectors V(:,i)'), pause
    clf
    I = 1;
    for i=[n-2,n-5,n-8,n-11]
      subplot(2,2,I); plot(1:n,XX(:,i)), axis([1,n,-1,1]);
      xlabel(['i = ',num2str(i)]), I = I + 1;
    end
    subplot(2,2,1)
    text(10,1.2,'Right generalized singular vectors XX(:,i)'), pause
    clf

    k_tsvd = 7; k_tgsvd = 6;
    X_I = tsvd(U,s,V,b,1:k_tsvd);
    X_L = tgsvd(UU,sm,XX,b,1:k_tgsvd);
    pause % Strike any key to continue
    subplot(2,1,1);
      plot(1:n,X_I,1:n,x,'x'), axis([1,n,0,1.2]), xlabel('L = I')
    subplot(2,1,2);
      plot(1:n,X_L,1:n,x,'x'), axis([1,n,0,1.2]), xlabel('L ~= I'), pause
    clf

    % Part 6.  No square integrable solution
    % --------------------------------------
    %
    % In the last example there is no square integrable solution to
    % the underlying integral equation (NB: no noise is added).
    %
    % Notice that the discrete Picard condition does not seem to
    % be satisfied, which indicates trouble!

    [A,b] = ursell(32); [U,s,V] = csvd(A);
    pause % Strike any key to continue
    picard(U,s,b); pause

    % This concludes the demo.
    echo off

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON