# Ill-conditioned Gram Matrix Assembly

I'm trying to find the best approximation to the function $e^x$ in the finite dimensional polynomial space $P_4$ with respect to the standard basis vectors $B=\{1,x,x^2,x^3,x^4\}$ with inner product $$(f,g)=\int_0^1 \frac{f(x)g(x)}{\left(x(1-x)\right)^{\frac{1}{2}}}dx$$

Let $w=\sum_{j=0}^4a_j\phi_j$ where $\phi_j=x^j$. Then, the best approximation of $w$ to $e^x$ satisfies the orthogonality property:

$$(e^x-w,\phi_i)=0$$ for $i=0,...,4$. Hence,

$$\sum_{j=0}^4 a_j (\phi_j,\phi_i) = (e^x,\phi_i)$$ for $i=0,...,4.$. This generates linear system of equations which I'm trying to solve. Unfortunately, the integrals are a bit difficult to evaluate by hand. So, I tried using a three point gaussian quadrature rule to approximate the integrals and solve the system. Unfortunately, I keep getting a nearly singular system. When I observed this, I tried using a higher order quadrature rule, but the matrix remains ill-conditioned.

I looked at my code over and over again and I can't seem to find any errors in the implementation. I'm tempted to believe that these inner products simply can't be evaluated well with a gaussian quadrature rule. Is this common for Gram Matrices? Should I use a different quadrature rule? Or is it a bug in my code?

I have provided my matlab code assembing this matrix system below for reference. Any help would be greatly appreciated.

n=5;
z1=-sqrt(3/5);
z2=0;
z3=sqrt(3/5);

x1=z1/2+1/2;
x2=z2/2+1/2;
x3=z3/2+1/2;

A=zeros(n,n);
b=zeros(n,1);

%assemble gramm matrix
for i=1:n
for j=1:n
F1=x1^(i-1+j-1)/sqrt(x1*(1-x1));
F2=x2^(i-1+j-1)/sqrt(x2*(1-x2));
F3=x3^(i-1+j-1)/sqrt(x3*(1-x3));
A(i,j)=(1/2)*((5/9)*F1+(8/9)*F2+(5/9)*F3);
end
G1=x1^(i-1)*exp(x1)/sqrt(x1*(1-x1));
G2=x2^(i-1)*exp(x2)/sqrt(x2*(1-x2));
G3=x3^(i-1)*exp(x3)/sqrt(x3*(1-x3));
b(i)=(1/2)*((5/9)*G1+(8/9)*G2+(5/9)*G3);
end


### Preliminary observations

• If you have a basis $\{\phi_i(x)\}$ of size $n$, you must have at least $n$ quadrature points. You can see this by writing the inner product as $(u,v) = u^T B^T W B v$ where $W$ is a diagonal weight matrix and $B_{ij} = \phi_j(x_i)$ is the basis evaluation matrix. If there are fewer than $n$ quadrature points, $B$ has fewer than $n$ rows, thus $\DeclareMathOperator{\rank}{rank} \rank(B^T W B) < n$.

• Your weight function cannot be integrated exactly using Gauss quadrature, depending on the order and desired accuracy, may have unacceptable errors. You can find papers on special quadratures for all sorts of specific functional forms. Note that there is also quadrature error integrating the right hand side. Fortunately, these functions are smooth and the weight effectively asks for more points near the endpoints, so Gauss quadrature is not that bad. See Trefethen (2008) Is Gauss Quadrature Better Than Clenshaw-Curtis for more on generic quadratures.

• The representation $A = B^T W B$ is a more convenient way to write your code.

• It is considered good Matlab style to use matrix and vector notation when convenient. It substantially compresses and simplifies the code.

### Refactored and working code

n = 5;
qpoints = 10;

% Quadrature points on [-1,1]; see Golub & Welsch (1969) and Trefethen (2008)
beta = .5 ./sqrt(1-(2*(1:qpoints-1)).^(-2)); % 3-term recurrence coefficients
T = diag(beta,1) + diag(beta,-1);            % Jacobi matrix
[V,D] = eig(T);                              % eigenvalue decomposition
[x,i] = sort(diag(D));                       % Gauss-Legendre points
w = 2*V'(i,1).^2;                            % weights
% Map quadrature to interval [0,1]
x = 0.5 + 0.5*x; w = 0.5*w;

% assemble Gram matrix
B = fliplr(vander(x,n));         % basis evaluation matrix
W = diag(w ./ sqrt(x.*(1-x)));   % diagonal quadrature weight matrix
A = B'*W*B;                      % Gram matrix
b = B'*W*exp(x);                 % exponential integrated against test functions

taylor = 1./factorial(0:n-1)';
ceoffs = [A \ b, taylor]


### Further observations

• The monomial basis is notoriously ill-conditioned. If you would like to work with higher degree polynomials, you should choose a well-conditioned basis such as (a) Legendre/Chebyshev polynomials, (b) Lagrange interpolants on good point sets like Gauss-Lobatto, or (c) Newton-form (divided differences) polynomials associated with a good point set. Of these, (b) is the most popular in PDE computations.
• The computed coefficients do not match the leading terms in the exponential because they are weighted-$L^2$ optimal on the interval $[0,1]$ rather than only in the limit as $x \to 0$. In contrast, the Taylor polynomials only apply to the $x\to 0$ limit and provide no bounds on a finite interval in any norm.
• You can produce optimal polynomials with respect to non-smooth norms such as $L^\infty(\Omega)$. Approximation in $L^\infty$, sometimes called "Chebyshev approximation", gives you optimal pointwise error bounds across the entire interval. This blog post has some nice figures, an introduction to the algorithms, and further links to the literature.
• 1st observation: I don't see how this is erroneous... $x^\alpha x^\beta\equiv x^{\alpha+\beta}$. What is it that I'm overlooking? 2nd: This must be the primary source of my ill-conditioned matrix, right? 3rd: Obviously, quadrature is inexact except for limited degree polynomials, but the error induced shouldn't lead to ill-conditioning in general, should it? 4th: Isn't it inefficient to formulate $B^TWB$ for higher order polynomials? 5th: The monomial basis always leads to a vandermonde matrix? 6th: I'm not sure how the weighted $L^2$ optimality affects the matching... could you explain? – Paul Oct 14 '12 at 3:53
• You're right, removed. I was thinking of them as general basis functions (which everything else works with). – Jed Brown Oct 14 '12 at 4:01
• 7th: What benefit can I get from using a non-smooth norm on this problem? – Paul Oct 14 '12 at 4:07
• 2. Monomials are the primary reason for ill-conditioning, provided the quadrature is sufficient. 3. Quadrature contributes some, compare cond(A) as you increase the quadrature order and if you replace with midpoint quadrature, for example. 4. The complexity is of the same order as your loops, but $B^TWB$ executes faster because it is vectorized. For specific integrals, you can use Clenshaw-Curtis which can integrate using FFT. 5. Yeah, but as mentioned before, usually you want to use a nicer basis. – Jed Brown Oct 14 '12 at 4:23
• 6/7. It depends what you want. $L^\infty$ provides pointwise error bounds on the interval. Point values can be arbitrarily bad in $L^2$ (depending on function regularity). If you want to approximate $\exp$ using a polynomial on an interval, you would very likely prefer to use $L^\infty$. This is sometimes called "Chebyshev approximation" in the literature. Note that Taylor polynomials do not usefully bound the error on any finite interval in any norm. – Jed Brown Oct 14 '12 at 4:27

The gram matrix of the powers is known to be very ill-conditioned. It is far better if you expand in terms of Chebyshev polynomials for an integral in [-1,1], and transform your integral to that range.