The answer is on the Wikipedia page on matrix mechanics, but I'll repeat some of it here. The intuition Heisenberg had was from classical mechanics, which is only valid in the correspondence limit, which is large values of the energy levels, and large n orbits.
There is a simple classical relation between the frequencies of the classical motion and the outgoing classical radiation (in the dipole approximation): the outgoing radiation intensity is in a periodic function, which produces outgoing periodic waves as the period of d(T), the dipole moment of the atom as a function of time. If you think the nucleus is stationary, this dipole moment is just the position of the electron, and the period is the orbital period.
The outgoing waves are just produced from the oscillating dipole moment, and they are long-wavelength compared to the atom, and they resolve into separated radiation at frequencies: 1/T, 2/T, 3/T, as any periodic function with period T does.
The classical function x(t) representing the position of a classical electron has a Fourier decomposition:
$$ x(t) = \sum_{k=-\infty}^\infty e^{ik\omega t} X_n $$
Where $X_n = X_{-n}^*$ (because X is real), and $\omega = {2\pi\over T}$. I am neglecting back reaction in this. The dipole moment is just the same thing times -e, and the field produced by this charge gives outgoing waves in a classical superposition of outgoing modes of the same frequencies $k\omega$, with an amplitude which goes as $X_n$.
If you have another classical function, say the momentum, you write it as a similar Fourier series:
$$ p(t) = \sum_k e^{ik\omega t} P_n$$
and the product $x(t) p(t)$ is, in Fourier space, the convolution of $X_n$ and $P_n$. This is what the "convolution algebra over the additive subgroup" means. The frequencies are integer spaced, with spacing equal to the inverse classical period and multiplying functions of time together is convolving their Fourier coefficients on this frequency set.
But in an atom, at large n,m, you know that the emitted photon energies are according to the formula:
$$ E_n - E_m = {C\over n^2 } - {C\over m^2}$$
Where C is some small multiple of the Rydberg constant. The frequency of the outgoing photon is just the energy (up to $\hbar$, which I take to be 1). These energies are not integer-spaced, so the frequencies, even at large n and m where you expect the motion to be nearly classical, are not exactly evenly spaced multiple of an inverse orbital period. This means that the classical orbital picture doesn't work, the frequencies are not classical.
If you assume n and m are big, but their difference is small compared to the absolute size of either:
$$ n = N$$
$$ m = N + k$$
Then the frequencies are approximately evenly spaced, since
$$ {1\over n^2} - {1\over m^2} \approx {-2k\over N^2} $$
And if you work out all the constants, this approximate even spacing is exactly the reciprocal period of the classical Bohr orbit corresponding to orbit N (or M, the difference is higher order). This was already known to Bohr, it is called the correspondence principle, and it was the major tool he used for finding the quantization rule. The correspondence principle says that the level spacing is the inverse of the classical period at large quantum numbers, and this allows you to derive the old quantum condition, that it is the action that should be quantized to leading semiclassical order.
Heisenberg is trying to give a dynamical theory which is quantum and complete, which contains a description like X(t) of the electron motion. He chooses to find the quantum analog of the time Fourier coefficients $X_k$. He knows that $|X_k|^2$ is the intensity of classical radiation emitted at frequency $k\omega$, and he also knows that the quantum mechanics thing is emitting with frequency, not $k\omega$, but $E_N - E_{N+k} \approx k\omega $$.
So he considers replacing the quantities $X_k$ with $X_{nm}$, according to the correspondence rule:
$$ X^N_k = X_{N (N+k}$$
which means, the classical Fourier series of the N-th Bohr orbit is the near-diagonal part of a matrix of quantities $X_{mn}$. The number of steps you are away from the diagonal is the order k.
The time evolution of Fourier coefficients (ignoring radiation) is simple:
$$ X_n(t) = e^{in\omega t} X_n(0)$$
From this, it is trivial to conclude that
$$ X_{nm}(t) = e^{i(E_n - E_m)t} X_n(0)$$
in quantum mechanics. This is the Heisenberg equation of motion.
The notion of product is convolution classically, so quantum mechanically, if you have $X_{nm}$ and $P_{rs}$, you want to produce a product of these which reduces to the classical convolution near the diagonal of the matrix. The products
$$ X_{nm} P_{rs} $$
have a frequency under time evolution which is $$ E_n - E_m + E_r - E_s $$, so if $m=r$, you get a cancellation and something which can be interpreted as a matrix element of a new quantity. Heisenberg immediately can guess:
$$ \sum_n X_{mn} P_{ns} = (XP)_{ms} $$
because this is what reduces to classical convolution in the correspondence limit. This is matrix multiplication, as Jordan and Born tell him. Next Heisenberg decides to find the form of the matrices from the known semiclassical Harmonic oscillator. This is very easy to guess, the classical motion has only one frequency, so the matrices are only 1 step off-diagonal.
Then Heisenberg translates the old quantum condition to matrix mechanics and finds the on-diagonal canonical commutation relation. With Born and Jordan, they find the full canonical commutation relation in the next paper. All this is from semiclassical reasoning, which is reproduced on Wikipedia in the way I found most clear.
Your specific concern
When discussing quantum systems, you construct the operator algebra from the known state-space using some operators which you already know what their state space is supposed to be. For a single particle, you start with the X operator basis where you have a continuous basis of x states. Then diagonalizing the Hamiltonian reproduces Heisenberg's picture.
If you are starting with 4 discrete states, you can make up any 4 by 4 matrices, and these will be some quantum observables. The correspondence between the index of the matrix and the spectral lines (differences of energy levels) is assuming that you are in a basis where the particle Hamiltonian is diagonalized --- that the states are energy states.
It is best if you read physics sources on this to learn it, not mathematics sources.