import numpy as np
import matplotlib.pyplot as plt


"""
Minimize f(x) = 0.5 x^T A x - b^T x
using steepest descent method with exact linesearch
"""

N = 20
# generate NxN random matrix
X = np.random.randn(N,N)
# generate NxN orthogonal matrix from it
Q = np.linalg.qr(X)[0]
# generate some random eigenvalues between lam_min and lam_max
lam_min = 1.0
lam_max = 300.0
lmbda   = lam_min + (lam_max-lam_min)*np.sort(np.random.rand(N))
lmbda[0]= lam_min
lmbda[-1]=lam_max
Lambda  = np.diag(lmbda)
A = np.matmul(Q,np.matmul(Lambda,Q.T))
# random vector
b = np.random.randn(N)
# A^{-1}b
xstar = np.linalg.solve(A,b)
# starting point
x0    = xstar + 1.0/lmbda[0]*Q[:,0] + 1.0/lmbda[-1]*Q[:,-1]


# function, derivative, analytical steplength
f = lambda x: 0.5*x.dot(A.dot(x)) - b.dot(x)
df= lambda x: A.dot(x) - b
steplength = lambda p: p.dot(p)/p.dot(A.dot(p))

# tolerance: stop when ||x-x_opt|| < tol
# in practice of course one stope when, e.g., |df(x)|<tol or similar
tol = 1.0E-06

x = x0
i = 0

err0 = np.linalg.norm(x-xstar,2)

while np.linalg.norm(x-xstar) > tol:
    # steepest descent direction
    p    = -df(x)
    err1 = np.linalg.norm(x-xstar,2)
    print("Iter: %3d, f=%15.6e, |grad f|=%15.6e, |x(k)-x*|/|x(k-1)-x*|=%15.6e" % \
          (i, f(x), np.linalg.norm(p,2), err1/err0))

    if np.linalg.norm(p) == 0.0:
        # success
        break
    # steplength: in this case we can determine it analytically
    alpha = steplength(p)
    x = x + alpha*p
    err0 = err1
    i = i+1

p    = -df(x)
err1 = np.linalg.norm(x-xstar,2)
print("Iter: %3d, f=%15.6e, |grad f|=%15.6e, |x(k)-x*|/|x(k-1)-x*|=%15.6e" % \
      (i, f(x), np.linalg.norm(p,2), err1/err0))

print("(lambda[n]-lambda[1])/(lambda[n]+lambda[1])) = %e" % ((lmbda[-1]-lmbda[0])/(lmbda[-1]+lmbda[0])))