Conjugate Gradient Optimizer

This routine lets you optimize large scale linear systems

% This example demonstrates the use of conjgrad.m
% The main advantage of conjgrad.m is that it takes handles to functions
% which perform the evaluation of the linear operator and its adjoint.
% The parameter space can be multidimensional.

%% Example #1 - matrix inversion

% generate well conditioned random matrix
N = 128;
[U,S,V]=svd(randn(N));
s=diag(S);
A=U*diag(s+max(s))*V;
b=randn(N,1);

% define the operator and its adjoint
operator = @(x) A*x;
adjoint = @(x) A'*x;
x0 = zeros(size(b));

res_limit = 1e-4;
max_steps = 100;
[x, Res] = conjgrad(x0, b, operator, adjoint, res_limit, max_steps);
plot(log10(Res));

%% Example #2 - deconvolution
N = 128;

% the convolution kernel is two dots. So the forward operator will make a
% twin-image shifted by 2 pixels.
kernel = zeros(N);
kernel(N/2,N/2) = 1;
kernel(N/2,N/2+2) = 1;
f_kernel = fft2(kernel);
test_image = randn(N);

% the adjoint of FFT is iFFT
% the adjoint of pointwise multiplication is multiplication by conjugate
operator = @(x) ifft2(f_kernel.*fft2(x));
adjoint = @(x) ifft2(conj(f_kernel).*fft2(x));
b = operator(test_image);
x0 = zeros(size(test_image));

res_limit = 1e-2;
max_steps = 100;
[x, Res] = conjgrad(x0, b, operator, adjoint, res_limit, max_steps);
plot(log10(Res))