Contents

% Solve the polynomial approximation problem with
% L-norm residual for L = 2, infty, and 1
clear all;

set(0, 'defaultaxesfontsize',18,'defaultaxeslinewidth',2,...
    'defaultlinelinewidth',2,'defaultpatchlinewidth',2)

generate data points

m = 20;
t = (-5:10/(m-1):5)';
b = 1./(t.^2+1);

add errors

b2 = b; b2(9) = b2(9) + 0.2;    % outlier noise
b3 = b.*(1 + 0.1*randn(m,1));   % 10% normally distributed noise
B = [b b2 b3];
Titles = ['noiseless case';
    'outliner noise';
    'Gaussian noise'];

Polynomial degree

n = 6;
A = ones(m,n);
for i = 2:n
    A(:,i) = A(:,i-1).*t;
end

loop through 3 noise sceniors

for j = 1:3

    b = B(:,j);
    disp(Titles(j,:));

noiseless case

outliner noise

Gaussian noise

L_2 approximation

    xl2 = A\b;
    fprintf('|| A x - b ||_2 = %12.6e \n\n',norm(A*xl2-b,2));

|| A x - b ||_2 = 5.896196e-01 

|| A x - b ||_2 = 6.201431e-01 

|| A x - b ||_2 = 6.450626e-01 

L_infty approximation

    AA   = [ A  -ones(m,1)
        -A  -ones(m,1)];
    cc   = [zeros(n,1); 1];
    bb   = [b; -b];

    [x,fval,exitflag] = linprog(cc,AA,bb);

    fprintf( 'linprog terminated with exitflag = %d \n', exitflag )
    fprintf( '        and objective function value   %12.6e \n', fval )

    xlinf = x(1:n);
    fprintf('|| A x - b ||_inf = %12.6e \n',norm(A*xlinf-b,inf));
    fprintf('eta               = %12.6e \n\n', x(n+1));

Optimization terminated.
linprog terminated with exitflag = 1 
        and objective function value   1.937218e-01 
|| A x - b ||_inf = 1.937220e-01 
eta               = 1.937218e-01 

Optimization terminated.
linprog terminated with exitflag = 1 
        and objective function value   1.937218e-01 
|| A x - b ||_inf = 1.937218e-01 
eta               = 1.937218e-01 

Optimization terminated.
linprog terminated with exitflag = 1 
        and objective function value   2.258630e-01 
|| A x - b ||_inf = 2.258630e-01 
eta               = 2.258630e-01 

L_1 approximation

    AA   = [ A  -eye(m)
        -A  -eye(m)];
    cc   = [zeros(n,1); ones(m,1)];
    bb   = [b; -b];

    [x,fval,exitflag] = linprog(cc,AA,bb);

    fprintf( 'linprog terminated with exitflag = %d \n', exitflag )
    fprintf( '        and objective function value   %12.6e \n', fval )
    x1 = x(1:n);
    fprintf('|| A x - b ||_1 = %12.6e \n',norm(A*x1-b,1));
    fprintf('sum y_i         = %12.6e \n\n', sum(x(n+1:end)));

Optimization terminated.
linprog terminated with exitflag = 1 
        and objective function value   1.946975e+00 
|| A x - b ||_1 = 1.946975e+00 
sum y_i         = 1.946975e+00 

Optimization terminated.
linprog terminated with exitflag = 1 
        and objective function value   2.146975e+00 
|| A x - b ||_1 = 2.146975e+00 
sum y_i         = 2.146975e+00 

Optimization terminated.
linprog terminated with exitflag = 1 
        and objective function value   2.128696e+00 
|| A x - b ||_1 = 2.128696e+00 
sum y_i         = 2.128696e+00 

plot results

    figure(j)
    plot(t,b,'bo'); hold on

    tt = -5:0.1:5;
    ff = 1./(tt.^2+1);
    plot(tt,ff,'k-'); hold on
    title(Titles(j,:));

    pp = xl2(n);
    for i = n-1:-1:1
        pp =  pp.*tt+xl2(i);
    end
    plot(tt,pp,'g--'); hold on

    pp = xlinf(n);
    for i = n-1:-1:1
        pp = pp.*tt+xlinf(i);
    end
    plot(tt,pp,'r:'); hold on

    pp = x1(n);
    for i = n-1:-1:1
        pp = pp.*tt+x1(i);
    end
    plot(tt,pp,'k-.'); hold off
    title(Titles(j,:));

    xlabel('x')
    legend('measurement','f','l_2 approx','l_\infty approx','l_1 approx')

loop end

end

Published with MATLAB® R2013b