function [xmin, ymin] = linesearch(Q, c, x0)
%LINESEARCH Steepest descent using line search.
%
%       1  T       T
%   y = - x Q x - c x
%       2

%   Ling Li
%   04/16/2001, Caltech

% make Q symmetric
Q = (Q + Q') / 2;

% The minimum should be
%       -1         1  T -1
%  x = Q  c, y = - - c Q  c
%                  2
xmin = inv(Q) * c;
ymin = - c' * xmin / 2;

x = x0; y = x'*Q*x/2 - c'*x;
xpath = [x]; ypath = [y];
while abs(y - ymin) >= 1e-7,
    % p is the steepest descent direction, i.e., - dy/dx
    p = c - Q*x;

    % a is the step length, decided by
    %             T
    %            p p
    %      a = -------
    %            T
    %           p Q p
    % since now p = - dy/dx.
    a = (p' * p) / (p' * Q * p);

    % iterate x: x_next = x + a p.
    x = x + a * p; y = x'*Q*x/2 - c'*x;
    % record the search path
    xpath = [xpath, x];
    ypath = [ypath, y];
end

% Calculate the condition number of Q
e = eig(Q);
k = max(e) / min(e);
% The convergence constant is [(k-1)/(k+1)]^2
C = (k-1)/(k+1); C = C * C;

% Let's estimate the convergence constant
dist = ypath - ymin;
%Ce = mean(dist(2:end) ./ dist(1:end-1));
% in fact, all entries in Ce is the same
Ce = dist(2) / dist(1);

disp(['condition number is ' num2str(k)]);
disp(['C = ' num2str(C) ', estimated = ' num2str(Ce)]);

% plot function and search path
[x1, x2] = meshgrid(xmin(1)+(-10:.1:10), xmin(2)+(-10:.1:10));
y = (Q(1,1)*x1.*x1 + Q(2,2)*x2.*x2)/2 + Q(1,2)*x1.*x2 - c(1)*x1 - c(2)*x2;
figure; hold on;
plot(xpath(1,:), xpath(2,:), 'o--');
for i = 1:size(xpath, 2),
    text(xpath(1, i)+0.2, xpath(2, i)+0.3, int2str(i-1));
end
contour(x1, x2, y, [ypath, ymin]);
title(['line search with k = ' num2str(k)]);
xlabel('x_1'); ylabel('x_2');

%figure; semilogy(ypath- ymin);
