# A Checking Gram-Schmidt

1. Oct 19, 2017

### joshmccraney

Hi PF!

I recently wrote a program in MATLAB to orthonormalize continuous functions on a given interval $x\in[a,b]$. Lets call the initial linearly independent functions I seek to orthonormalize $L_i$ and let's call these functions $\psi_i$ after I orthonormalize them.

After running my script, I perform an integration check, where I create an $N\times N$ matrix $A:A_{ij} = \int_a^b \psi_i\psi_j\,dx$, $N$ being the number of $L_i$'s I use. Theoretically $A$ should be the identity matrix.

When I let $N=5$, my maximum error outside the diagonals for a particular sequence of functions (too difficult to explain which functions) is very small $O(10^{-5})$. However, when I let $N=20$, the maximum error off the diagonals is $O(1)$.

My question is, do any of you know a set of functions that, once orthonormalized, give another analytic set of functions? I thought about taking $L_i = x^i:x\in[-1,1]$ to generate the Legendre polynomials, but these polynomials are not orthonormal on $x\in[-1,1]$. I used $\sin(nx)$ and $\cos(nx)$ and they work, but not all sequences of functions seem to (like the first function I eluded to in the previous paragraph).

Does anyone know either what the Legendre polynomials are orthonormal over OR a set of linearly independent functions that have a known analytic corresponding orthonormal set of functions?

I'm trying to check my code and this would be very helpful. For your info, here's the script I wrote
Code (Text):

N = 20;% number of functions

dx = 0.001;
x = -pi:dx:pi;
Li = zeros(N,size(x,2));
for n = 1:size(Li,1)
Li(n,:) = cos(n*x)./sqrt(pi);% linearly independent functions to orthonormalize
end

psi = zeros(size(Li));
ortho = zeros(size(psi));

% gram-schmidt
for ii = 1:size(psi,1)
for s = 1:ii-1
ortho(ii,:) = ortho(ii,:)+psi(s,:)*trapz(psi(s,:).*Li(ii,:))*dx;
end% for s
psi(ii,:) = (Li(ii,:)-ortho(ii,:))./sqrt(trapz((Li(ii,:)-ortho(ii,:)).^2)*dx);% Gram-Schmidt on Li to construct psi
end% for ii

% check orthogonality

for kk = 1:size(Li,1)
for jj = 1:size(Li,1)
orBC(jj,kk) = trapz(psi(kk,:).*psi(jj,:))*dx;% orthogonality check on psi (should be Identity matrix)
end% for jj
end% for kk

for kk = 1:size(Li,1)
for jj = 1:size(Li,1)
orBCexact(jj,kk) = trapz(Li(kk,:).*Li(jj,:))*dx;% orthogonality check on Li (should be Identity matrix)
end% for jj
end% for kk

2. Oct 20, 2017

### NFuller

I'm not completely sure what you are asking but the Legendre polynomials are orthogonal over the interval $[-1,1]$ such that
$$\int_{-1}^{1}P_{n}(x)P_{m}(x)dx=\frac{2}{2n+1}\delta_{nm}$$

3. Oct 20, 2017

### joshmccraney

Thanks for replying! Yea, I saw that orthonormality relationship. So after closer inspection I think the script I attached works for most functions (it certainly does for Legendre polynomials. However, I have an interesting set of functions, which are linear combinations of Legendre polynomials over $[-1,1]$, which, when I enter them in as $L_i$, they are not all orthogonal.

Let me explain my problem in detail and if you don't mind you can help me out. Let's call the $ith$ Legendre polynomial $p_i:i=0,1,2...$. I have a piecewise function $f$ defined as
$$f = \begin{cases} f_1 = \sum_{i=0}^N b_i p_i(x)&x\in [-1,\zeta_1] \\ f_3 = 0 & x \in (\zeta_1,\zeta_2)\\ f_2 = \sum_{i=0}^N c_i p_i(x)&x\in [\zeta_2,1] \end{cases}$$
I have three boundary conditions $f$ is subject to. In matrix form, these three boundary conditions form a $3\times2(N+1)$ matrix $M$ that operates on the vector coefficients $\mathbf{b} = [b_0,b_1...b_N]^T$ and similarly for $\mathbf{c}$. In matrix form this looks like $$[M]\begin{bmatrix} \mathbf{b}\\ \mathbf{c} \end{bmatrix} = \mathbf{0}.$$
To solve this I compute $null(M)$, which is a $2(N+1)\times 2N-1$ matrix. Then there are $2N-1$ coefficient vectors $[\mathbf{b}; \mathbf{c}]^T$ that satisfy the boundary condition matrix $M$. So, once computing $null(M)$, I know $\mathbf{b}, \mathbf{c}$, and thus I know $f_1$ and $f_2$ and consequently I know $f$. Since there are $2N-1$ coefficient vectors, I actually have $2N-1$ linearly independent $f$ functions. However, these functions $f$ are not necessarily orthonormal over $[-1,1]$. This is where I apply Gram-Schmidt. When I apply Gram-Schmidt, it will transform $f_i\to\psi_i:i=0,1,2...2N-1$ where $\psi_i$ should be orthonormal.

To verify each $\psi_i$ is othonormal over $x\in[-1,1]$, I make a matrix $A : A_{ij} = \int_{-1}^1 \psi_i \psi_j\,dx$. If everything is correct, $A = I_{2N-1}$, the identity matrix. When $N = 5$ everything works out very well, and $A \approx I$. However, when $N=10$, $A \neq I$, where early entries, say $A_{3,4}$ or $A_{2,7}$ are approximately zero yet higher entries $A_{17,18}\approx 1$. I should say $A_{ii}=1$ for any $N$ selected.

When troubleshooting this, I was curious about the vector coefficients from $null(M)$. So I plotted $\psi_{15-19}$ and noticed they are almost identical. Any idea what's going on here? Attached are three plots: $\psi_i:i=0,1,...,19$, $\psi_i:i=15,...,19$, and a zoomed in view of $\psi_i:i=15,...,19$, all plotted over $x$.

Sorry, I keep forgetting how to post a picture directly on this post.

#### Attached Files:

• ###### psi.pdf
File size:
361.3 KB
Views:
34
4. Oct 20, 2017

### jasonRF

EDIT: I was looking at your code again, and it looks like may be already using the modified Gram-Schmidt? If so then my first guess is wrong. But to be honest I am confusing myself trying to look at your code - a problem with my brain at the end of a long day, not your code!

One issue that might be causing you grief is that, even for the "simple" case of finite vectors (not functions) the classical Gram-Schmidt algorithm is numerically unstable on a floating point machine. You are essentially implementing classical Gram-Schmidt so that could be the problem; your use of trapz for estimation of the inner products may also be adding to the instability? Anyway, see the Wikipedia page:

https://en.wikipedia.org/wiki/Gram–Schmidt_process

Your symptom of having higher order $\psi$ functions being very similar is the same thing I have seen before when Gram-Schmidt is unstable. You might would not see the instability if you start with a nice set of functions on a nice interval (say, sinusoids on a [-pi,pi] or Legendre polynomials on [-1,1]) except for very large N, but for some choices of functions you will see it for N relatively small, like 10.

jason

Last edited: Oct 20, 2017
5. Oct 20, 2017

### joshmccraney

This is exactly the advice I was looking for. I'll do some research, thanks!

6. Oct 23, 2017

### joshmccraney

So I used the algorithm attached on wikipedia, and for most functions I get the same result as my code in post 1. However, when I apply this to the function I refer to in post 1, (which I call $f$ in the code below) the plot looks nasty for large $\psi$ functions. The largest value of $A$ off the diagonal is 0.19, so that is better, but still not small enough to consider it orthogonal. Attached is a picture of the new $\psi$ values for the function.

Any idea what's happening? My code in post 1 seems to be more stable, as the "orthonormal" functions appear to be continuous, yet obviously my code gives linearly dependent solutions for large $N$. What's going on here? I'm happy to provide any additional info that would help.

I could program this in Mathematica and avoid all numerical errors, since the function $f$ that's giving me trouble is a linear combination of Bessell functions, but it would be so nice to get this working in MATLAB. Plus, this is a good learning experience for me.

Code (Text):

V = f';
n = size(V,1);
k = size(V,2);
U = zeros(n,k);
U(:,1) = V(:,1)/sqrt(trapz(V(:,1).*V(:,1))*dx);
for i = 2:k
U(:,i) = V(:,i);
for j = 1:i-1
U(:,i) = U(:,i) - ( U(:,i)'*U(:,j) )/( U(:,j)'*U(:,j) )*U(:,j);
end
U(:,i) = U(:,i)/sqrt(trapz(U(:,i).*U(:,i))*dx);
end
U = U';

#### Attached Files:

• ###### psi.pdf
File size:
58.7 KB
Views:
31
7. Oct 23, 2017

### jasonRF

One thing about that code is that you are not consistent about when you use a trapz inner product (defining U(:,1) and normalizing U(:,i) at the end of each for i= ... loop) versus a 'normal' inner product (inside the for j= ... loop).

Have you tried running this using simple inner products first instead of your trapz inner products?

jason

8. Oct 23, 2017

### joshmccraney

Shoot, good call. So I've fixed the Gram-Schmidt as shown below but it's still giving me very discontinuous outputs nearly identical to those shown in post 6. I did run this on simpler Legendre polynomials and it works out well. Really confused why it's not working out here. I could send you a plot of the functions that are not orthonormalizing under this procedure if you'd like?

Code (Text):

V = f';
n = size(V,1);
k = size(V,2);
U = zeros(n,k);
U(:,1) = V(:,1)/sqrt(trapz(V(:,1).*V(:,1))*dx);
for i = 2:k
U(:,i) = V(:,i);
for j = 1:i-1
U(:,i) = U(:,i) - (trapz(U(:,i).*U(:,j))*dx )/( trapz(U(:,j).*U(:,j))*dx )*U(:,j);
end
U(:,i) = U(:,i)/sqrt(trapz(U(:,i).*U(:,i))*dx);
end
U = U';

9. Oct 24, 2017

### jasonRF

I'm looking at your code and that is definitely not the modified Gram-Schmidt - it is the unstable classical Gram-Schmidt. My fault for not making sure wikipedia was correct!

The stable modified Gram-Schmidt (here coded up as a QR decomposition) would look more like,
Code (Text):
% QR decomposition via modified Gram-Schmidt.  V = Q*R
[m,n] = size(V);  % columns of V have initial vectors
R = zeros(n,n);   % R is upper triangular matrix
Q = zeros(m,n);   % columns of Q will have orthonormal vectors
for k=1:n
R(k,k) = norm(V(:,k));
Q(:,k) = V(:,k)/R(k,k);
for jj=k+1:n
R(k,jj) = Q(:,k)'*V(:,jj);
V(:,jj) = V(:,jj) - Q(:,k)*R(k,jj);
end
end
A trapz version would only have to modify two lines.

In any case, I recommend looking at your vectors with standard tools and standard (not trapz) inner products to make sure your problem is well conditioned. One way is to use Matlab's built-in QR decomposition. I suspect it does not use Gram-Schmidt in any form, but instead uses either Householder reflectors or Givens rotations, which are even better than the modified G-S, especially if orthogonality of your vectors is of prime concern. For your problem try,
Code (Text):
% let columns of V be vectors to make orthonormal
[Q,R] = qr(V,0);  % columns of Q give orthonormalized set of vectors

If you try this and the matrix Q is orthogonal (with respect to standard inner products) to an acceptable degree, then try to use the modified Gram-Schmidt above and see how it does. If it is fine, then a trapz version of the modified Gram-Schmidt code should be adequate for you. If not, then you will need to decide if the trapz stuff is really that much better than a standard inner product. If it is, then you will need to code up something like a Householder algorithm with trapz, which will require you to re-derive parts because all of the versions you are likely to find will implicitly be assuming standard inner products. This is do-able, but will take some real work on your part.

Hope that helps.

Jason

EDIT: any book on numerical linear algebra should cover all of this. Examples are Golub and Van Loan, "Matrix Computations", and Watkins "Fundamentals of Matrix Computations". Any edition for either book.

Last edited: Oct 24, 2017