Archives

Archives

April 10 Be sure to read over the solutions to A4 and A5 in preparation for the test. Receiving a perfect score on a problem does >>not<< mean that you've nailed all the underlying issues.

April 5. We will use Wednesday's class (4/9) as a review session

April 4. Prelim on Friday, April 11. Topics = everything since prelim 1. I.e., Assignment 4 (eigenvalues) and Assignment 5 (n=1 nonlinear root finding and optimization.).

April 2 CS 422 Final : Monday, May 12, 9-11:30AM, Hollister 314. Let me know if you have an exact conflict or 3 exams in a 24hr period including ours. We cannot give early exams because of airplane tickets, early summer jobs etc.

April 2. Test scripts for A5 are up.

Mar 26 Polygon Smoothing scripts/solutions/insights here.

Mar 25 Things to do for P3. (a) When you apply power method with M and a start vector with centroid zero, how close is the kth iterate to the space spanned by the c and s vectors? Relevant: find scalars alpha and beta to minimize

norm( [c s]*[alpha;beta] - x^{k})

where x^{k} is the k-th iterate, normalize so its 2-norm is one.

(b) Run ShowP with the initial x and y vectors as linear combos of c and s.

(d) Spot patterns in M^k*c/norm(M^k*c and M^k*s/norm(M^k*s)

Mar 20 Although it turns out not to make a difference, the original x and y vectors in ShowP are not normalized. So, for the sake of clarity the line

x = rand(n,1); y = rand(n,1);

should be changed to

  x = rand(n,1); y = rand(n,1); x = x - sum(x)/n; y = y -sum(y)/n;

The on-line version of ShowP has been updated accordingly.

Mar 11 Some comments on the Q's that are generated for each Jacobi update. The c-s computation for the 2-by-2 case is discussed in the text. It involves solving a quadratic which of course has two roots. You should choose that root which corresponds to the smaller of the two possible rotation angles. This means that the [c s;-s c] that you generate should be as close as possible to eye(2,2). This turns out to be important for convergence. It guarantees that the 2-by-2 rotations are converging to eye(2,2). Otherwise there are examples where large off-diagonal elements are "chased around" rather than zeroed. Take this into account when you design ParJac.

For the BlockJac it is more complicated because now a 2m-by-2m orthogonal matrix Qij is associated with each Jacobi subproblem. If you compute Qij using the Matlab schur function you will observe that a large number of sweeps are required. This is because the Qij's are not converging to eye(2m,2m). To rectify this, you should try to make Qij look as much as possible as eye(2m,2m) . One way to ``clean up'' the Qij produced by schur is to replace it by Qij*P where P is a permutation that moves the largest row entries to the diagonal position . If these maximal entries are negative, multiply the corresponding column by -1.

Q = 
   -0.2507    0.4556   -0.8542
   -0.9653   -0.0514    0.2559
    0.0726    0.8887    0.4527

you should rearrange it to

Q = 
  0.8542   0.2507    0.4556   
 -0.2559   0.9653   -0.0514    
 -0.4527  -0.0726    0.8887       

to make it look more like the identity,  This maneuver will still give 
you a Q that diagonalizes the subproblem.
    
Mar 11  Test scripts 
for A4 are up.  
Mar 7  Assignment 4 
has been handed out.
Mar 7  Prelim 1 solution 
available. median = 85. 


 
 
Feb 28 I added some hints to recent probOfDays
Feb 28 Office hours today 10:30-Noon. You can pick-up 
your A3. After noon, Cindy in 4146 Upson has them.
Feb 28 A3's are graded. Median = 16. Main mistakes:
P1; You have to ordert he unknowns righth if you are to live off 
of one QR factorization: a1, b1, a2, b2, etc. If you dont do this you cannot 
simply grab Q(:,1:2*n) and R(1:2n,1:2n) and expect to get the right coefficient 
values.
P2. Sqrt(theta) is required!
P3. When confronted with an underdetermined system Cx = d, where 
C is m-by n with m<n, you cannot assume that C(1:m,1:m) is nonsingular. I.e., x 
= [C(1:m,1:m)\d ; zeros(n-m,1)] may not exist.
P4. Minimizing norm( [x-x0 y-y0]*c - (z-z0)) and then setting  
a = [c(1) c(2) -1] is not the same optimization problem.
Feb 27 A3 solutions available.    
Feb 25 More Hints
P1. If you have the thin QR factorization of an m-by-n matrix A 
you also have the thin QR factorization of A(:,1:r) where r<n.
P2. Givens QR for an m-by-n A goes like this
  for j=1:n
     for i=m-1:-1:j
         Zero A(i+1,j) with a Givens rotation and update A
     end
  end

Also, you do not need to explicitly form the Q-matrix to solve 
the LS problem.
P3. Chapter 2-type of hint. To solve a square linear system 
M1*M2*x = b where M1 and M2 are nonsingular, you first solve M1*y = b for y and 
then solve M2*x = y for x. Of course, you need factorizations of M1 and M2 to 
pull it off.
Feb 25 The test script P1 needs an m =1000 at the start. 
The online version has been corrected.
Feb 21 The test scripts for A3 are now
available. Some hints and comments:
P1. Do not compute any unnecessary factorizations. Also, just 
submit the numerical output. The plots are just for you to observe the quality 
of the fit for various choices of n.
P2. You are not allowed to use the backslash operator except to 
solve square, triangular linear systems, Apply the Givens rotations efficiently. 
P3.  The problem should have stated "make effective use of 
the Matlab qr function." And BTW, you'll need to call qr twice. There is a typo, 
the solution vector x is in R^{r} not R^{n}.   
P4. Express the problem in matrix-vector terms.
>> Points will be deducted for poorly commented code. In 
particular, write good specifications for the four required functions
Feb 20 Assignment A3 is posted.  
Test scripts will be available on Friday. The dues date/time is "hard" because 
solutions have to be posted immediately to give you time to studey for the 
Friday (2/29) quiz. More on the Quiz later.
Feb 18  The explanation for why (just about) 
everybody is getting extremely high efficiency measures in P4  remains a 
mystery! I will try to figure it out by Wednesday. One theory was that backslash 
is not "smart" if handed a triangular matrix. I tried this benchmarking script 
on my computer and it confirms that if you double n, work in A\b goes up roughly 
by 4 if A is triangular. What does it say about your computer/Matlab version?
clc
nRepeat = 300
disp(' ')
disp('Time A\b when A is triangular')
disp(' ')
for n = [100 200 400 800]
    A = triu(rand(n,n));
    b = rand(n,1);
    tic
    for k=1:nRepeat
        x = A\b;
    end
    t = toc;
    fprintf('%4d  %10.6f\n',n,t) 
end

disp(' ')
disp('Time A\b when A is full')
disp(' ')
for n = [100 200 400 800]
    A = rand(n,n);
    b = rand(n,1);
    tic
    for k=1:nRepeat
        x = A\b;
    end
    t = toc;
    fprintf('%4d  %10.6f\n',n,t) 
end


 
 
Feb 14   Either re-download load P3 or change 
the lines
      v = d;

      u = y*sigma;
to
      u = d;

      v = y*sigma;
You'll get better cosines this way!
Feb 12 There is a missing "=" in the equation at the top 
of the second page. It should be in between the two vectors.
Feb 12   There remains some "off by 1" confusion in P1 despite 
my correction  below. To clear this up, here is a complete specification 
for BandLU:
  function [c,d,f] = BandLU(a)
% a is a column n-1 vector
% c is a column n-1 vector
% d is a column n vector
% f is a column n-1 vector
% If A = eye(n,n) + diag(a,1) - diag(a,-1), L = eye(n,n) + diag(c,-1)
% and U = eye(n,n) + diag(f,1), then A = LU
 
Also, there was a mistake in P1. Download it again

 
Feb 11 Assignment 2 is posted. Couple of typos in the 
handout associated with P1:

	There is a stray "c"  in the displayed n=5 L-factor
	The (3,3) entry in the displayed n= 5 U-factor should be d_3
	In the prose just beneath that equation, should have 
	a(1:n-1), not a(1:n).

Feb 4. Couple of things about A1. In P2, please print out 
eps*cond(A) not eps*relativeError. The idea is for you to see correlations 
between the relative error and this product. Note that you should be able to 
deduce the exact solution from how the system is set up. In P3 your solution 
will involve matrix-matrix products where one of t he matrices is a permutation. 
You may regard these as zero-flop operations. (In practice, it amounts to an O(n^2) 
data motion calculation.)
Jan 28 First Assignment is 
available.
Jan 20. A reasonable "first approximation" of the course 
is the fall 
2006 addition of CS 421. There are some syllabus changes however. In 
particular, we will NOT be doing the initial value problem or any other 
differential equation problems. The fall course CS 421 (a.k.a. Math 421) covers 
that material. Note that the two courses in this year-long sequence may be taken 
in either order.