% A comparision of gradient descent and conjugate gradient (guess who wins) function main() % data A=[17, 2; 2, 7]; % the matrix b=[2, 2]'; % right-hand side x0=[0, 0]'; % the initial guess % linewidth and font size lw= 2; fs = 25; % colors red=[0.867 0.06 0.14]; blue = [0, 129, 205]/256; green = [0, 200, 70]/256; black = [0, 0, 0]; white = 0.99*[1, 1, 1]; % Set up the plotting window figure(1); clf; set(gca, 'fontsize', fs); hold on; axis equal; axis off; s = 0.16; x = A\b; Ax = x(1)-s; Bx = x(1)+s; Ay = x(2)-2.0*s; By = x(2)+s; plot([Ax Bx Bx Ax Ax], [Ay Ay By By Ay], 'color', blue, 'linewidth', lw/2); % plot a blue box s=0.005; plot(Ax-s, Ay-s, '*', 'color', white); plot(Bx+0.5*s, By+0.5*s, '*', 'color', white); %markers Box = [Ax Bx Ay By]; axis (Box); % plot the contours of the quadratic form associated with A and b plot_contours(A, b, Box, lw, blue); % Do conjugate gradient and gradient descent. % For the first one, start a bit shifted so that the two graphs don't overlap. shift = 0.0015*[1, -1]; small_rad=0.002; tol = eps; x = conj_gradient(A, b, x0, tol, lw, red, small_rad, shift); x = grad_descent (A, b, x0, tol, lw, green, small_rad); % text small = 0.015; text(x0(1)-2*small, x0(2)-1.6*small, 'x', 'fontsize', fs); text(x0(1)-0.5*small, x0(2)-3*small, '0', 'fontsize', floor(0.7*fs)); text(x(1)+small, x(2)+small, 'x', 'fontsize', fs); % some balls for beauty small_rad = 0.003; ball(x0(1)+shift(1)/2, x0(2)+shift(2)/2, small_rad, blue) ball(x(1), x(2), small_rad, blue) % save to disk as eps and svg saveas(gcf, 'Conjugate_gradient_illustration.eps', 'psc2'); plot2svg('Conjugate_gradient_illustration.svg'); function x = conj_gradient(A, b, x, tol, lw, color, small_rad, shift) r=A*x - b; d=-r; while norm(r) > tol % a pretty ball for beauty, to cover imperfections when two segments are joined ball(x(1)+shift(1), x(2)+shift(2), small_rad, color); alpha = -dot(r, d)/dot(A*d, d); x0 = x; x = x + alpha*d; r=A*x - b; beta = dot(A*r, d)/dot(A*d, d); d0 = d; d = -r + beta*d; plot([x0(1), x(1)]+shift(1), [x0(2), x(2)]+shift(2), 'color', color, 'linewidth', lw) end function x = grad_descent(A, b, x, tol, lw, color, small_rad) r=A*x - b; d=-r; while norm(r) > tol % a pretty ball for beauty, to cover imperfections when two segments are joined ball(x(1), x(2), small_rad, color); alpha = -dot(r, d)/dot(A*d, d); x0 = x; x = x + alpha*d; r=A*x - b; beta = 0; %beta = dot(A*r, d)/dot(A*d, d); d0 = d; d = -r + beta*d; plot([x0(1), x(1)], [x0(2), x(2)], 'color', color, 'linewidth', lw) end function plot_contours (A, b, Box, lw, color); N=200; % number of points (don't make it big, code will be slow) E = A\b; % the exact solution, around which we will draw the contours B = 0.12; [X, Y]=meshgrid(linspace(Box(1)-B, Box(2)+B, N), linspace(Box(3)-B, Box(4)+B, N)); % X and Y coordinates % the quadratic form f= (1/2)*x'*A*X-b'*x; f = inline('0.5*A(1, 1)*X.*X + A(1, 2)*X.*Y+0.5*A(2, 2)*Y.*Y-b(1)*X-b(2)*Y', 'X', 'Y', 'A', 'b'); Z = 0.5*A(1, 1)*X.*X + A(1, 2)*X.*Y+0.5*A(2, 2)*Y.*Y-b(1)*X-b(2)*Y; % prepare to draw the contours x0 = A\b; f0 = f(x0(1), x0(2), A, b); No = 25; % number of contours Levels = (linspace(f0, 1, No)-f0).^2+f0; % Plot the contours with 'contour' in figure(2), and then with 'plot' in figure(1). % This is to avoid a bug in plot2svg, it can't save output of 'contour'. figure(2); clf; hold on; for i=1:length(Levels) figure(2); [c, stuff] = contour(X, Y, Z, [Levels(i), Levels(i)]); [m, n]=size(c); if m > 1 & n > 0 % extract the contour from the contour matrix and plot in figure(1) l=c(2, 1); x=c(1,2:(l+1)); y=c(2,2:(l+1)); figure(1); plot(x, y, 'color', color, 'linewidth', lw/2); end end figure(1); function ball(x, y, r, color) Theta=0:0.1:2*pi; X=r*cos(Theta)+x; Y=r*sin(Theta)+y; H=fill(X, Y, color); set(H, 'EdgeColor', 'none');