Inverse of an infinitely large matrix?

Question

This is probably a trivial problem for some people, but I've spent quite some time on it: What is the inverse of the infinite matrix $$ \left[\begin{matrix} 0^0 & 0^1 & 0^2 & 0^3 & \ldots\\ 1^0 & 1^1 & 1^2 & 1^3 & \ldots\\ 2^0 & 2^1 & 2^2 & 2^3 & \ldots\\ 3^0 & 3^1 & 3^2 & 3^3 & \ldots\\ \vdots & \vdots & \vdots & \vdots &\ddots \end{matrix}\right] $$ (Assume that $0^0=1$ for this problem).

I'm not sure if this problem has a solution or is well-defined, but if it has a solution, it would help greatly in a ton of stuff I'm doing (mostly related to generating functions and polynomial approximations). I began by taking the inverse of progressively larger square matricies, but I didn't see any clear pattern.

Answer 1

The linear transformations of $\Bbb R^\infty$ Take the form of row-finite matrices (all rows only have finitely many nonzero entries) or column finite matrices (analogously) depending on if you're writing the vectors on the right or left of transformations.

This matrix is neither row nor column finite, so it can't represent a linear transformation, invertible or not.

Answer 2

Numerically, we get some very interesting results for the matrix $M_{ij}=(i-1)^{j-1}$ if it is expressed in the Fourier basis with alternating row signs, $$A := \begin{bmatrix}1 \\ & -1 \\ && 1 \\&&& -1 \\ &&&& \ddots \end{bmatrix}\mathcal{F}^{-1} \left[\begin{matrix} 0^0 & 0^1 & 0^2 & 0^3 & \ldots\\ 1^0 & 1^1 & 1^2 & 1^3 & \ldots\\ 2^0 & 2^1 & 2^2 & 2^3 & \ldots\\ 3^0 & 3^1 & 3^2 & 3^3 & \ldots\\ \vdots & \vdots & \vdots & \vdots &\ddots \end{matrix}\right] \mathcal{F}.$$

As the discretization size of $M$ goes up, $M$ becomes extremely ill-conditioned and the entries blow up, so numerical calculations using standard doubles (15 digits of accuracy) will fail when it becomes much larger than 10-by-10. However, there are toolboxes that let you do computation with greater precision, and I used one with 512 digits of accuracy to produce the following images of the real and complex entries of the matrices $A$ up to size 128-by-128:

The color of the $(i,j)$'th pixel in each plot represents $A_{ij}$, the value of the $(i,j)$'th entry of $A$. The first row is the real part of $A$, and the second row is the imaginary part. Red means large positive value, blue means large negative value, and the maximum real value in the plot is shown in the middle. The picture is big but scaled down for display on math.stackexchange - you can open it in a new tab to see it in more detail.

It looks like in this Fourier basis the normalized matrices are converging to an integral operator with a smooth kernel, $$\frac{1}{N}A v \rightarrow \int_0^1 (K(x,y) + iJ(x,y)) v(y) dy,$$

where $K$ and $J$ are the smooth functions in the pictures, and $N$ is some renormalization factor.

Since the kernel is smooth, the action of $A$ will annihilate highly oscillatory functions. Recalling the definition $A = \text{diag}(1,-1,\dots) \mathcal{F}^{-1} M \mathcal{F}$, we see exactly how $M$ is ill-conditioned and what functions are in its numerical null space - functions that are the Fourier transform of a highly oscillatory functions.

Turning this around, there should exist an inverse $A^{-1}$ for the renormalized limit, acting on a space of functions that are sufficiently smooth. Indeed, the following is a plot of the spectrum of the renormalized $A$ in the 128-by-128 case (top right), and it's dominant singular vectors (left, real on top imaginary on bottom):

---

Here's the Matlab code I used:

%Using Advanpix multiprecision computing toolbox, http://www.advanpix.com/
mp.Digits(512+9);
mp.GuardDigits(9);

jjmin = 2;
jjmax = 7;
jjrange = jjmax - jjmin + 1;

for jj=jjmin:jjmax
    N = 2^jj;

    %Generate original matrix M, where M_nm = (n-1)^(m-1)
    v = mp((0:N-1)'); 
    M = mp(zeros(N,N)); 
    for kk=0:(N-1) 
        M(:,kk+1)=v.^kk; 
    end

    %Generate matrix D F^(-1) M F, where F is the fft, and 
    %D is the diagonal matrix with diagonal [1, -1, 1, -1, ...]
    FMF = mp(zeros(N,N)); 
    for k=1:N 
        ek = mp(zeros(N,1)); 
        ek(k)=1; 
        FMF(:,k) = ifft(ifftshift(M*fft(fftshift(ek)))); 
    end
    for k=1:N 
        FMF(k,:) = (-1)^(k-1)*FMF(k,:); 
    end

    %plot it
    subplot(2, jjrange,jj - jjmin+1)
    imagesc(real(double(FMF)))
    title(['N=', num2str(N)]);

    subplot(2, jjrange, jjrange + jj -jjmin+1)
    imagesc(imag(double(FMF)))
    format short
    title(['max= ',num2str(max(real(FMF(:))),3)])

end

subplot(2,jjrange,1)
ylabel('real(D F^{-1} M F)')

subplot(2,jjrange,jjrange + 1)
ylabel('imag(D F^{-1} M F)')