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: # 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])))