Re: 'orth' command -- question

[John B. Thoo]

Hi, James.

On Apr 17, 2012, at 5:25 AM, James Sherman Jr. wrote:

> On Tue, Apr 17, 2012 at 8:05 AM, John B. Thoo <address@hidden> wrote:
> Hello, again, James.
> 
> On Apr 16, 2012, at 11:21 PM, John B. Thoo wrote:
> 
> > Hello, James.
> >
> > On Apr 16, 2012, at 11:04 PM, James Sherman Jr. wrote:
> >
> >> On Tue, Apr 17, 2012 at 1:37 AM, John B. Thoo <address@hidden> wrote:
> >> Hi.  I'm trying to understand the command 'orth'.
> >>
> >> Example 1
> >> ---------
> >> octave-3.2.3:46> A = [1, 2, 2; 2, 1, 2; 2, 2, 1];
> >>
> >> octave-3.2.3:47> [V, LAMBDA] = eig (A); P = orth (V)
> >> P =
> >>
> >>  0.62060  -0.53058   0.57735
> >>  0.14920   0.80275   0.57735
> >> -0.76980  -0.27217   0.57735
> >>
> >> octave-3.2.3:48> P'*A*P
> >> ans =
> >>
> >> -1.0000e+00   2.7756e-17  -9.4369e-16
> >>  1.1102e-16  -1.0000e+00   8.6042e-16
> >> -8.8818e-16   7.7716e-16   5.0000e+00
> >>
> >>
> >> So, it appears that 'orth' provides an orthonormal basis of eigenvectors 
> >> of  A.
> >>
> >> Example 2
> >> ---------
> >> octave-3.2.3:66> A = [4, 1, 0; 1, 4, 1; 0, 1, 4];
> >> octave-3.2.3:67> [V, LAMBDA] = eig (A); P = orth (V)
> >> P =
> >>
> >> -0.023793   0.865699   0.500000
> >> -0.588348   0.392232  -0.707107
> >> -0.808257  -0.310998   0.500000
> >>
> >> octave-3.2.3:68> P'*A*P
> >> ans =
> >>
> >>  4.9791e+00  -6.5271e-01   2.2204e-16
> >> -6.5271e-01   4.4351e+00   9.9920e-16
> >>  4.4409e-16   6.1062e-16   2.5858e+00
> >>
> >> Now it appears that 'orth' does _not_ provide an orthonormal basis of 
> >> eigenvectors of  A.
> >>
> >> Why does 'orth' appear to behave differently in the two examples?
> >>
> >> Thanks.
> >>
> >> ---John.
> >>
> >> -----------------------------------------------------------------------
> >> "Ten thousand difficulties do not make one doubt....  A man may be annoyed 
> >> that he cannot work out a mathematical problem ... without doubting that 
> >> it admits an answer."
> >>
> >> ---John Henry Newman [_Apologia_, p. 239 in Project Gutenberg's
> >>  <http://www.gutenberg.org/ebooks/22088>]
> >>
> >>
> >> _______________________________________________
> >> Help-octave mailing list
> >> address@hidden
> >> https://mailman.cae.wisc.edu/listinfo/help-octave
> >>
> >> Hi John,
> >>
> >> I'm not quite sure what you mean by an "orthonormal basis of eigenvectors 
> >> of A".  Specifically in your second example, the matrix V is already 
> >> orthonormal (V'*V = eye(3) and its columns are composed of eigenvectors of 
> >> A), so I'm not sure what your intention is in calling orth is for.  Also, 
> >> since orth just looks for an orthonormal basis for the column space of V, 
> >> and since the eigenvectors found have distinct eigenvalues, they span all 
> >> of R^3, thus any orthonormal basis for R^3 will suffice.  So, I'd say that 
> >> it was just luck involved that the first case worked as you had expected 
> >> and the other case didn't, and it just depends on the particular numerical 
> >> algorithm involved.
> >
> > Thanks for your reply.
> >
> > My question was really whether one could use 'orth' "blindly," that is, 
> > without having to check first if 'eig' already returned an orthonormal set 
> > of eigenvectors.  If I understand correctly, one cannot avoid having to 
> > think. :-)  Also, now I see that 'orth' returns an orthonormal basis that 
> > has the same span only, but does not intentionally respect eigenspaces. 
> > (Right?)
> >
> > The reason for seeking an orthonormal set of eigenvectors of  A is to 
> > orthogonally diagonalize  A.  I'm teaching an introductory linear algebra 
> > course and, for the first time, I'm assigning some problems that require 
> > the use of numerical computation software.  (I figured I needed to creep 
> > into the 21st century finally.)  The textbook's exercises refer to the 
> > "Dark Side," but I've asked my students to use Octave.  My problem is that 
> > I'm still very much a novice at using Octave, and am trying to stay 1/2 
> > step ahead of my students. :-O  So far it's been great fun!
> 
> Forgetting about eigenvectors, if given an orthogonal matrix  A,  apparently 
> 'orth (A)' does not necessarily return  A.  Why?  Thanks.
> 
> ---John.
> 
> -----------------------------------------------------------------------
> "Ten thousand difficulties do not make one doubt....  A man may be annoyed 
> that he cannot work out a mathematical problem ... without doubting that it 
> admits an answer."
> 
> ---John Henry Newman [_Apologia_, p. 239 in Project Gutenberg's
>   <http://www.gutenberg.org/ebooks/22088>]
> 
> 
> 
> I certainly understand where you're coming from.  I teach a linear algebra 
> course as well (and it just so happens we're getting to diagonalization 
> tomorrow).  At least the book I use tends to just refer to "numerical 
> computation software" and not the Dark Side. :)

The textbook we're using (I won't tell you that it's Kolman/Hill) has a section 
in the back that introduces the Dark Side specifically.


> You actually run into this behavior a lot with software, because the idea of 
> "=" is very flaky in a mathematical sense when talking about floating point 
> numbers, and there are lot of tricks involved to make the output numerically 
> stable for large matrices.  

That's very interesting.


> I don't know the specific algorithm behind orth, so I can't give any 
> particular insight into this algorithm.  Another example that I've found that 
> has this type of behavior is LU factorization (which drives my students 
> nuts), because even for small matrices, Octave will perform row switching 
> before attempting to factorize into L and U matricies (again for numerical 
> stability purposes), so that the L and U that Octave gives satisfies:
> LU = PA
> where P is a permutation matrix, and not the factorization they were 
> assigned, LU = A.

Yes, I noticed that, too, with the LU factorization earlier in the semester.


> I know this isn't exact answer to your question, but I hope it helps.

Yes, it helps.  I'm not only still very much a novice at using Octave, but I'm 
quite inexperienced with using numerical computation software of any kind.  The 
types of issues you mentioned (like numerical stability) are not ones that come 
to mind immediately.  Now I know that I need to be careful and that things may 
not turn out as I expect.  It's a good lesson, and one for me to share with my 
students.

Thanks for your help and patience.

---John.

-----------------------------------------------------------------------
"Ten thousand difficulties do not make one doubt....  A man may be annoyed that 
he cannot work out a mathematical problem ... without doubting that it admits 
an answer."

---John Henry Newman [_Apologia_, p. 239 in Project Gutenberg's 
   <http://www.gutenberg.org/ebooks/22088>]