CAAM 454/554 Homework Assignments

The problems are numbered as in the textbooks. TB refers to Trefethen and Bau, and NW to Nocedal and Wright (2nd Ed.).

Problems starting with (554) will be required only from CAAM 554 students (similarly for those with (454)).

HW 10: Due 4/19, Tue. (554 Solution)

(40 points) Assignment: Primal-Dual Interior-point method for LP
Submit your output (truncated the iteration history if too long) along with your code and a short write-up (typed) explaining what you observe.
(554, 5 points each) Problems NW 13.5, 14.8, 14.9, 14.11. (Correction: in problem 14.9, the equation (14.10) should be (14.16))

HW 9: Due 4/7, Thur. (554 Solution)

(15 points) Write a matlab function for BFGS method with the following interface:
```
function [x,sflag,hist] = myBFGS(func,x0,tol,maxit)
%
% Quasi-Newton method using inverse BFGS update 
% (if you do backtracking, make sure you check y'*s > 0)
%
% Usage: 
%   [x,sflag,hist] = myBFGS(func,x0,tol,maxit);
%
% INPUT:
%       func - matlab function for f and g evaluations
%         x0 - starting point (column vector)
%        tol - stopping tolerance for gradient norm
%      maxit - maximum number of iterations
%
% OUTPUT:
%          x - computed approximate solution
%      sflag - status flag, contains number of iterations
%              if convergence achieved or 0 if maxiter reached
%      hist(1,:) - vector of function values at all iterations
%      hist(2,:) - vector of gradient norms  at all iterations
```
--- Use the identity matrix for the initial inverse Hessian approximation.
--- Stop when the norm of gradient is less than the tolerance.
--- Download the instructor's code zBFGS.p (also pre-2007 version zBFGS07.p), the generalized Rodenbrock function (see page 98, Ex. 4.3) rosenbrock_g.p, and the test script test9.m, run it for n = 10, 100, 200.
--- Make your output format similar with that of the instructor's code.
--- Download the script test9b.m, run it with your codes myNewton, myBFGS, myCDG2 on the function fhw8.m, all using the stopping criterion that the gradient norm is less than the tolerance. See a diary of output from a run by the instructor diary9b.txt.
--- Submit the output (truncated printout and graphs) along with your code and a short write-up (typed) explaining what you observe.
(15 points) Implement the Gauss-Newton method with a line search to solve a set of nonlinear least squares problems. The interface of the code should be:
```
 
    [x,iter] = MyGN(func,x,tol,maxiter,prt,varargin);
```
where func is a function name that returns [r,J], prt is a switch for printing iteration info or not, and varargin is a variable list of input arguments to be used by function func. The stopping criterion should be norm(g) < tol, where g = J'*r is the gradient vector.
Download and unzip the set of codes in nls.zip. Download the instructor's codes ZyGN.p and ZyLM.p (and pre-07 versions ZyGN07.p and ZyLM07.p). Study and run the test code testNLS.m. Without input, it runs the instructor's codes (CodeRange = 3:6) by default. We included 33 test problems, as described in the paper MGH.pdf. The test problems have already been coded. A program FdNewton.m is included that implements Newton's method for solving nonlinear systems using a finite-difference scheme to approximate Jacobian. You may start your programming by studying and modifying this program.
Submit your code and printout from the test program testNLS.m for the following two runs:
```
        testNLS(1, 1.e-4,1:4);   (easy run) 
        testNLS(10,1.e-8,1:4);   (hard run)
```
(554, 20 points) Problems NW 6.1(a), 6.3, 6.4, 6.9 (prove Psi(B)>= n instead).
(554, 15 points) Problems NW 10.1, 10.4, 10.6.

HW 8: Due 3/30, Wed. 6pm (slide under DH3090 door) (Solution (a)) (Solution (b))

(10 points) Write a gradient (steepest) descent Matlab function:
```
	function [x,iter,nf] = myGD(func,x0,tol,maxit,varargin)
```
to minimize a function defined by the function handle func starting from initial guess x0. You interface func by
```
	[f,g] = feval(func,x,varargin{:}); % function and gradient
	f = feval(func,x,varargin{:});     % function value only
```
Use back-tracking line search and choose your own parameters. The stopping criterion is either
```
	crit := norm(g)/(1 + abs(f)) < tol;
```
or the maximum iteration number maxit is reached. The outputs include an approximate solution, the iteration number taken, and the number of function evaluations taken. You should mimic the printout format as in the instructor's code.
Download the test scripts testDenoise.m and testDeblur.m, the functions TVL2N.m and TVL2B.m, and the instructor's code zGD.p (also for pre-2007 versions zGD07.p). You can run the test scripts without your codes in place.
Run the test scripts with your code (you will need image processing toolbox). Submit your code and printout including graphs, and a short report (typed) explaining what you observe in these experiments.
(10 points) Modify your GD code to have a function that can also use nonlinear CG methods:
```
	function [x,iter,nf] = myCGD(func,method,x0,tol,maxit,varargin)
```
where method = 1,2,3,4 corresponding to GD, CG-PR, CG-PR+ and CG-FR, respectively. Check the search direction at each iteration. If a direction is not descent, use either the opposite direction or the steepest descent direction instead. Download the test scripts testDenoiseCGD.m, testDeblurCGD.m and the instructor's code zCGD.p. You should mimic the printout format as in the instructor's code.
Run the scripts with your code. Submit your code and the printout including graphs, and a short write-up (typed) explaining what you observe in these experiments.
(10 points) Write a Matlab function to minimize C² convex functions using Newton's method (without line search):
```
	function [x,iter,nf] = myNewton(func,x0,tol,maxit,varargin).
```
Stop the iterations when the norm of the gradient is less than the tolerance. Test the solver on f(x) = (x'*A*x)^2/4 + (sum(x) - n)^2/20. You need to supply a Matlab function [f,g,H] = fhw8(x,A) to calculates the function, gradient and Hessian values for this function.
Change the stopping criterion in myCGD.m so that it stops when the relative change in function value becomes less than the tolerance. Call it myCGD2.m. Download the test script testNewton.m, and run it for n = 500, 1000. Submit a short write-up (typed) explaining your observations in these experiments, such as convergence rates, with printout and other necessary supporting evidence.
(554, 10 points) Study Theorems 3.5, 11.2 and 11.3, understand what they say, their proofs and connections between them. Then stete and prove a theorem, following the model of Theorem 11.3, for the damped Newton method that solves the nonlinear equations r(x)=0. Specifically, (11.17) is replaced by
```
	x_k+1 = x_k + (1 - η_k)p_k, 	for some η_k in [0,η],
```
where p_k is the Newton direction at x_k.
(554, 10 points) Consider the following algorithm for solving nonlinear systems. At the current iterate x_k, do a Newton iteration to generate an intermediate point. Then do another Newton-like iteration from this intermediate point but using the same Jacobian (evaluated at x_k) to get the next iterate. (a) Write the iteration into a fixed-point iteration form x_k+1 = S(x_k). (b) What is the rate of convergence for this algorithm (under suitable conditions)? Sketch a formal argument for your claim. (This part will not be graded. Just do your best.)

HW 7: 3/17 Thur. (total points: 454/554 = 40/60) (Solution)

All functions are assumed to be from Rⁿ to R unless otherwise specified.

Show that the intersection of a collection of convex sets is convex.
Show that if f(x) is a convex function, then the set {x: f(x) &le 0} is a convex set.
Prove that the point-wise maximum of a (finite or infinite) set of convex functions is convex.
Prove that the set of real symmetric positive semidefinite matrices is convex.
(554) Prove that the maximum eigenvalue is a convex function over the set of real symmetric positive semidefinite matrices.
(554) Prove that a C² function is convex if and only if its Hessian is positive semidefinite at all x.

HW 6: 3/10 Thur. (Solution) (Solution2)

Sec. 5.6, page 15: Exercises 1, 2, 3 in the Krylov-Subspace Methods Notes. (554) Exercise 4.
Write a Matlab function function X = solve6(A,B,tol,maxit) to solve linear systems in a matrix form AX + XA = B, where all involved are square matrices. Assume that A is symmetric positive definite, and use the Matlab pcg function to solve the system by supplying a function (handle) to perform matrix-vector multiplications corresponding to the linear operation AX + XA. Turn off the on-screen printout of pcg function by specifying 2 outputs (see help pcg). You may use Matlab function reshape to convert matrix to vector and vice versa.
Download the test script test_hw6.m and the instructor's code zsolve6.p. You can run the test script without or with your code solve6.m. First try on small n values (less than 100), then try on a large value (500 or larger). Submit your code, and output (on screen and graph, if any) from a small and a large n value. Your code should have a reasonable speed compared to the instructor's code. Write a paragraph to describe your observations.

(454) Problems NW (2nd ed.) 2.1, 2.2, 2.9.
(554) Problems NW (2nd ed.) 2.13 - 2.15.

HW 5: Due 2/22, Tue. (total points: 454/554 = 30/40) (Solution TB38.2) (Solution TB38.6)

Program a Matlab function: function [x,iter] = myCGN(A,b,tol) that applies CG to the normal equations of Ax = b where A is nonsymmetric. Download the test script testCGN.m and the function mat4gmres.m. Run it for a large n, and submit your code and printout, and a short write-up (typed) explaining (not recounting) what you observe in this experiment.
Problems TB 38.2, 38.6.
(554) Exercises 1-3 (page 12) in the Krylov-Subspace Methods Notes.

HW 4: Due 2/10, Thur. (total points: 454/554 = 20/30) (Solution)

Modify your (non-restarting) GMRES code into function [x,iter] = myGMRES1p(A,b,tol,maxit,M1,M2) so that it utilizes pre-conditioners M1 and M1 whenever they appear in the input list. The pre-conditioner is the product M = M1*M2 and the pre-conditioned system is inv(M)Ax = inv(M)b. See Matlab "help gmres" for more clues.
Download the test script testGM1p.m and the instructor's code zGMRES1p.p. You can run the script without your code. When your code is ready, run the test script for n = 500 or larger (but use smaller values in the testing phase). Submit a copy of your program, the printput (graph not necessary), and a short write-up (typed) explaining what you observe in this experiment.
Modify your (non-restarting) GMRES code into a restarting GMRES code function [x,iter] = myGMRES2(A,b,p,tol,maxit) so that it restarts every p iterations. The input maxit is the maximum number of outer iterations allowed, and the output iter has 2 components, iter(1) is the number of outer iterations taken and iter(2) the number of inner ones taken in the last outer iteration. See Matlab "help gmres" for more clues.
Download the test scripts testGM2.m and testGM2a.m, and the instructor's code zGMRES2.p. You can run these scripts without your code. When your code is ready, run testGM2 for m = 100 or larger (as large as you see fit), and also run testGM2a. Submit a copy of your program, the printputs (graph included) from your runs, and a short write-up (typed) explaining what you observe in the experiment (testGM2a).
Let A = aI + UV^T be non-singular and diagonalizable where U and V are m by r with r < m, and a is a scalar. For this kind of matrices A, prove that GMRES converges in at most r+1 iterations.
(554) Let A = [e₂ e₃ ... e_m e₁] and b = e₁, where e_j is the j-th column of the identity (so A is a permutation of the identity). (a) Prove that GMRES requires m iterations to solve Ax = b by looking at the Krylov subspaces generated by (A, b). (b) Derive the characteristic polynomial of A and then its eigenvalues.

HW 3: Due 2/3, Thur. (total points: 454/554 = 30/40) (Solution)

Implement the (non-restarting) GMRES in Matlab function [x,iter] = myGMRES1(A,b,tol,maxit) where the output iter is the iteration number taken. Use the zero initial guess and the same stopping criterion as before. You should solve the inner least squares problems in an efficient way (but start your implementation with a simple one-liner and replace it later). The implementation given in the instructor's notes is recommended.
Download the test script testGM1.m, the function fdm4pde1.m for matrix generation, and the instructor's GMRES code zGMRES1.p. You can run the test script with or without your code. After you have written your code, start running the test script for a small m (say, m = 40) to debug and improve your code.
Submit your code and all printout for runs with m = 60, 80, or with m = 80, 100 (or even larger) if you can handle them comfortably.
Problems TB 33.1, 35.5.
(554) Problems TB 33.2(b)-(e).

HW 2: Due 1/27, Thur. (total points: 454/554 = 20/40) (Solution)

function [x,iter,info] = mySOR(A,b,omega,maxit,tol)

iter

omega

Download the test script testsor.m and the instructor's code zSOR.p. You can run "testsor.m" with or without your code.
After you have debugged your code, run "testsor.m" with mode = 1 and N <= 100 (try smaller N first), and then run it with mode = 2 and N >= 200 (or larger), and with a good omega values that you can estimate. Submit a copy of your program, the screen printout and the pictures for the above 2 runs.

Prove in detail that GS converges for any initial guess if the matrix A is symmetric positive definite.
(554) Excercises 2 and 4 in the lecture notes on Stationary Iterative Methods.

HW 1: Due 1/20, Tue. (total points: 454/554 = 30/45) (Solution)

Write two Matlab functions:

function [x,info,res] = myJacobi(A,b,maxit,tol) function [x,info,res] = myGS(A,b,maxit,tol)
to implement the Jacobi and Gauss-Seidel methods, respectively, where maxit is the maximum number of iterations. Use the zero vector to start. The stopping criterion is that the residual 2-norm, norm(b-A*x), becomes less than tol*(1 + norm(b)). The output info is a character string indicating the status of termination: either "Converged at iteration ?" or "Maxit reached". And res is a vector containing residual norms, norm(b-A*x), for all iterations, counting iterations from 1.
Download the instructor's codes zJacobi.p, zGS.p and the test script test1.m. Read the test script.
Run test1.m for N = 80, 120, 160 (or larger if so desired), and with both method = 'Jacobi' and method = 'GS'. Submit your programs and output from these runs.
Prove in detail that the Jacobi method converges for any initial guess if the matrix A is column strictly diagonally dominant.
(554) Prove that the scheme x <-- x + a*(b-Ax) converges for any initial guess if the matrix A is symmetric positive definite and a>0 is less than 2 divided by the maximum eigenvalue of A.
(554) Prove that in any matrix norm ||.||, ||A^k||^1/k ---> &rho(A) as k goes to infinity.