Least Squares and Fourier Series

Question

I have a little bit of problem figuring out the relation between Fourier series and Least Squares.

As far as I understand, LS is a way of minimizing the quadratic error between a measured value $y_i$ and a linear combination of functions evaluated at points $x_i$. So, in order to solve for least squares we have to have the $(x_i,y_i)$ and the functions $f_j(\cdot)$ that are going to be linearly combined.

If I would like to approximate a set of points $(x_i,y_i)$ using Fourier basis functions, I could assemble a matrix as

$$ \begin{pmatrix} 1 & \cos(\omega x_1) & \sin(\omega x_1) \\ 1 & \cos(\omega x_2) & \sin(\omega x_2) \\ \vdots & \vdots & \vdots\\ 1 & \cos(\omega x_n) & \sin(\omega x_n) \\ \end{pmatrix} \begin{pmatrix} A_0 \\ A_1 \\ A_2 \\ \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{pmatrix}$$

To obtain the values of the coefficients $A_0,A_1,A_2$ using least-squares I can solve for $x = (A^TA)^{-1}A^Tb$

Now, my doubt arises since it seems to me that the coefficients $A_0,A_1,A_2$ computed with the least-squares formula can be found with the Fourier transform. Are the coefficients $A_0,A_1,A_2$ the same as the Fourier coefficients $C_k$ ?

EDIT

clc; close all;

N = 10;
x = 1:N;
%y = rand(1,N);
y = [0.7099 0.6745 0.7374 0.1144 0.7732 0.0021 0.0357 0.4804 0.8107 0.0410];

% Polynomial Fit
A = [(x.*x)' x' ones(N,1)];
c = A\y'

% Fourier Fit - K+1 Coefficients
A = [];
K = 8;
for k = 0:K
    A = [A exp(-1j*2*pi*k*x/N)'];
end
c = A\y'
c = (fft(y)/length(y))'
% Test the approximation on a bigger range
x_nodes = linspace(0,10,1000);
f = c(1);
for i = 2:K
    f = f + c(i)*exp(1j*2*pi*i*x_nodes/N);
end

If I try to use the complex form of the Fourier Series $ f = C_0 + \sum_k C_ke^{-j2\pi kx/L}$ and follow the least-squares approach and compare the results to the coefficients of the FFT, the coefficients seem to be different. What am I missing?

Answer 1

In a word, yes.

A least squares problem arises in the context of an over determined system. That means that you have more equations than unknowns and would be the case if you had more data points than Fourier coefficients that you're trying to find. This would give you the "best fit" Fourier series of a given order. For example, you could find the best fit of a 4 term Fourier series to a set of 20 data points.

On the other hand, the discrete Fourier transform of a set of points always gives the same number of Fourier coefficients as input points. Thus, if you start with 20 points you will get 20 Fourier coefficients. Additionally, the first N Fourier coefficients are exactly the same as a least squares fit of a Fourier series with only N terms.

The key here is that the Fourier basis is an orthogonal basis on a given interval. The math works out so that the least squares best fit based of a lower order Fourier series is exactly equivalent to the truncated FFT. See this PDF for more details.

To verify, you can easily try this yourself. Pick a set of data points and compute the DFT and the least squares best fit and compare them. Here's an example in MATLAB:

    clear

%interval length and number of points
L = 10;
N = 15;

%data points
x = (linspace(0,L,N+1))'; x(end)=[];
y = x.^2; %function values

% FFT:
yk = fft(y)/length(x); %divide by length(x) for normalization

% pick out the first half of the coefficients
% (coefficients are symmetric for real valued functions)
yk = yk(1:ceil((length(x)+1)/2));

% Transform coefficients to the form:
%       y = A0 + sum( Ak*cos(...) + Bk*sin(...) )
A0 = yk(1); %constant term
A = 2*real(yk(2:end)); %cosine coefficients
B = -2*imag(yk(2:end)); %sine coefficients

%account for difference in highest wave number for odd/even numbers of points
if mod(numel(x),2)==0
    A(end) = A(end)/2;
end

%reconstruct the original function for verification
yy=A0;
for k=1:length(A)
    yy = yy + A(k)*cos(2*pi/L*k*x) + B(k)*sin(2*pi/L*k*x);
end



%least squares fit, unknowns C = [A0 A(1) B(1) A(2) B(2) ... ]'

%number of sin() and cos() coefficients to use:
N = 4;

%set up matrix:
M = zeros(length(x),1+2*N);
M(:,1) = 1;
for k=1:N
    M(:,2*k) = cos(2*pi/L*k*x);
    M(:,2*k+1) = sin(2*pi/L*k*x);
end

C = M\y; %equivalent to (M'*M)\(M'*y)

% least squares coefficients:
A0L2 = C(1);
AL2 = C(2:2:end);
BL2 = C(3:2:end);

%reconstruct the function values
yyL2=A0L2;
for k=1:length(AL2)
    yyL2 = yyL2 + AL2(k)*cos(2*pi/L*k*x) + BL2(k)*sin(2*pi/L*k*x);
end

% Plot:

figure(1)
clf
hold on
plot(x,y,'ko')
plot(x,yy,'b-')
plot(x,yyL2,'r--')
legend('data','Fourier Series','Least Squares','location','northwest')
title('Comparison of function values')
xlabel('x')
ylabel('y')

figure(2)
clf
hold on
plot(A,'bo')
plot(B,'bs')
plot(AL2,'r.')
plot(BL2,'rx')
plot(0,A0,'bo')
plot(0,A0L2,'r.')
legend('Fourier Series, A_k','Fourier Series, B_k','Least Squares, A_k','Least Squares, B_k','location','north')
title('Comparison of coefficient''s values')
xlabel('wave number, k')
ylabel('coefficients')

This code produces the following plots. The second plot in particular shows that the coefficients of a truncated Fourier series are the same as the least squares best fit.

Answer 2

They're related, but not the same, for a square system, in which the least squares fit becomes an interpolation. You have the relationship that

$c_k = \pm (a_k+ib_k)/2$ if $k\neq 0$

and $c_0 = a_0$ for $k=0$. You can do this by noting that the Fourier expansion over a period of $\tau$ gives

$f(t) = a_0 + \sum_{k=1}^\infty a_k \cos(\frac{2\pi k}{\tau}t) + \sum_{k=1}^\infty b_k \sin(\frac{2\pi k}{\tau}t)$.

Rewriting using Euler's formula gives

$f(t) = \sum_{k=1}^\infty c_k e^{i(2\pi k/\tau)t}$

with the above relationship between $c$ and $a,b$ above.

I'm not sure if the analogy holds for least squares.