Why does NumPy's irfft2 of rfft2 lead to a matrix with one less column when the original matrix has an odd second index?

In order to make sure that irfft2 is in fact the inverse of rfft2, you need to let it know the exact shape of your input data when reversing the transformation.

Like so:

import numpy as np
x = np.ones((4, 3))
ix = np.fft.rfft2(x)
rx = np.fft.irfft2(ix, x.shape)

This is necessary precisely for the reason you highlight with your question: The way the transformed data (the "spectrum", ix in your example) is represented for real-valued input data (x) depends on whether the number of samples is odd or even in any of the dimensions.

The (i)rfft* family of functions are all tailored to the common use case where the input data is a series of real numbers, i.e. not complex numbers. The discrete Fourier transform of such an input is usually complex-valued, but has a special symmetry: the negative-frequency components are the complex conjugates of the corresponding positive-frequency components. That is, the spectrum contains essentially the same numbers twice, and half the spectrum already contains the information necessary to reconstruct the input data. Which makes sense: The spectrum is a series of complex numbers, which can be represented as two real numbers each, but the input data does not have that "complexity", as it is real-valued.

Then again, "half the spectrum" is not that clear a term when the length of the data (and thus of the full spectrum) may be odd or even. Mathematically, these two cases must be treated slightly differently. Which is why the length of the data is needed when reconstructing the input signal.

As the NumPy documentation of rfft notes for the one-dimensional case:

If n is even, [the last array element of the spectrum] contains the term representing both positive and negative Nyquist frequency (+fs/2 and -fs/2), and must also be purely real. If n is odd, there is no term at fs/2; [the last array element of the spectrum] contains the largest positive frequency (fs/2*(n-1)/n), and is complex in the general case.

And the documentation of irfft further explains:

The correct interpretation of the hermitian input depends on the length of the original data, as given by n. This is because each input shape could correspond to either an odd or even length signal. By default, irfft assumes an even output length which puts the last entry at the Nyquist frequency; aliasing with its symmetric counterpart. By Hermitian symmetry, the value is thus treated as purely real. To avoid losing information, the correct length of the real input must be given.

So an even-length signal is the default. Which is why you only run into this issue for odd lengths of the array dimension. The documentation of irfftn notes specifically that it is the inverse of rfftn only if called like irfftn(rfftn(x), x.shape).

FFTs in general do not handle odd-length inputs. They actually want powers of 2. When you FFT an odd-length vector of reals, you lose some information. If you try your experiment with a (4,4), you'll see that the output exactly matches the input.