Convergence of the regularized form of the Column-Kaczmarz algorithm.

Demonstrates the change in the rate of convergence of the regularized form of the Column-Kaczmarz algorithm.


Author & Reference

Ivanov Andrey, Graduate student,

Dept. of Applied Mathematics, S. P. Korolyov Samara State Aerospace University (National Research University), Faculty of Computer Science, Samara, Russia.

See more for Regularization Tools dependency

addpath ('../libs/regu/');

See more for AIRTools dependency

addpath ('../libs/AIRtools/AIRtools/');

Generate test problem

Discretization of the "famous" Fredholm integral equation of the frst kind deviced by D. L. Phillips, see more in [6].

[A,f,u_true] = phillips(512); %or may be shaw(512);
[m,n] = size(A);

Add noize to f vector

noize_value = 0.015;
noize = randn(size(f));
noize = (noize/norm(noize)) * norm(f) * noize_value;
d = (norm(noize) / norm(f)) * 100;

fn = f + noize;

N = [1:1:6000];

Execute the Kaczmarz alg

[x0 info0] = kaczmarz(A,fn,N);
error_x0 = zeros(1,size(x0,2));
for i = 1:1:length(error_x0)
    error_x0(i) = norm(x0(:,i)-x0(:,end),2);

Some estimation perturbation

delta = norm(f-fn);

One of the heuristic methods to select regularization parameter

See more: V. A. Morozov, Methods of solving incorrectly posed problems, Springer Verlag, New York, 1984.

sigmaMax = max(svd(A));
alpha = ((sigmaMax^2)*delta/(norm(f)+delta));

options.stoprule.type = 'NR';
options.stoprule.epsilon = 10^-8;
options.alpha = alpha;

Execute the Kaczmarz-Column-Regularization (KCR) alg.

[x1 info1] = columnkaczmarz(A,fn,N,[],options);

error_x1 = zeros(1,size(x1,2));
for i = 1:1:length(error_x1)
    error_x1(i) = norm(x1(:,i)-x1(:,end),2);

Execute the Kaczmarz-Random-Column-Regularization (KRCR) alg.

[x2 info2] = randcolumnkaczmarz(A,fn,N,[],options);

error_x2 = zeros(1,size(x2,2));
for i = 1:1:length(error_x2)
    error_x2(i) = norm(x2(:,i)-x2(:,end),2);

Execute the QO-Kaczmarz-Column-Regularization (QO-KCR) alg.

[x3 info3] = qoptimalcolumnkaczmarz(A,fn,N,[],options);

error_x3 = zeros(1,size(x3,2));
for i = 1:1:length(error_x3)
    error_x3(i) = norm(x3(:,i)-x3(:,end),2);

Plot the results

loglog (1:1:length(error_x0), error_x0, 'b',...
          1:1:length(error_x1), error_x1, 'r',...
          1:1:length(error_x2), error_x2, 'k',...
          1:1:length(error_x3), error_x3, 'y');
legend('KM', ...
       'KCR', ...
       'KRCR', ...

plot(1:1:length(x0(:,end)), x0(:,end), 'b', ...
     1:1:length(x1(:,end)), x1(:,end), 'r', ...
     1:1:length(x2(:,end)), x2(:,end), 'k', ...
     1:1:length(x3(:,end)), x3(:,end), 'y', ...
     1:1:length(u_true)   , u_true   , 'c');
legend('KM', ...
       'KCR', ...
       'KRCR', ...
       'QO-KCR', ...
       'Exact Solution');