Conjugate Gradient Optimizer

This routine lets you optimize large scale linear systems

Description

% 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))

• MATLAB release: MATLAB 7.9 (2009b)

[Bruno Luong]

Not sure what advantage of conjgrad over Matlab builtin function LSQR (where preconditioning is possible):

x = lsqr(@afun,b);
    
function y = afun(x,transp_flag)
        if strcmp(transp_flag,'transp')
            y = A'*x;
        else
            y = A*x;
        end
    end

[Peter]

Well, the only advantage is that the vectors x and y don't necessarily have to be ordered in a row format, they can be multidimensional arrays. The internal function mydot can handle their dot products.

Otherwise you are right, there is no much difference.

The pcg matlab function is also an alternative.