My textbook writes that
The net energy stabilization due to the occupation of a bonding MO is equal to the net energy destabilization due to the occupation of the corresponding antibonding MO.
Essentially, it claims that these two quantities are equal.
(Image source: SparkNotes)
Shouldn't the antibonding orbital be destabilized more than the bonding orbital is stabilized, though?
TL;DR
"The destabilisation of an antibonding MO is equal to the stabilisation of a bonding MO, relative to the constituent AOs."
Several key assumptions are inherent in this statement.
If we make all five assumptions, then the statement above is correct.
If we only make assumptions 1 through 4, then MOs will still have the same form as above (i.e. everything else about the MO diagram is correct), but the statement above will be incorrect: the antibonding orbital will be more destabilised than the bonding orbital is stabilised.
The question linked in the comments does give a qualitative explanation, so depending on what you are looking for, it may be worth checking out.
The MOs that you see are one-electron wavefunctions that are obtained by looking at the single-electron species $\ce{H2+}$. Finding the exact stationary states of this system (by analytical solution of the Schrödinger equation) is a bit complicated, so the strategy is to find a set of approximate stationary states for the electron in $\ce{H2+}$, which we can label as $\{|i\rangle\} (i = 1, 2, 3 \cdots)$. These states are called molecular orbitals.
For multi-electron molecules, like $\ce{H2}$, a further approximation is usually invoked, in expressing the electronic wavefunction of $\ce{H2}$ as a properly antisymmetrised product of these states $\{|i\rangle\}$. Therefore, the MOs, which were only "rigorously" derived for $\ce{H2+}$, may be extended to molecules such as $\ce{H2}$. Very loosely speaking, we may assign one stationary state $|i\rangle$ to one electron, which corresponds to placing an electron in an orbital. Because of the Pauli exclusion principle, each stationary state may only be assigned to a maximum of two electrons.
These are obviously crude approximations and it is a surprise that it even works decently well. In particular the issue of electron-electron repulsions in multi-electron systems is glossed over; in this very simple model, the MOs, being derived from a one-electron system, are entirely unaffected by them, which doesn't really make sense. As the comments on the question have alluded to, there are much more sophisticated ways of obtaining MOs nowadays, but since the claim in the question is derived using this model, we have to use it.
Once we have assumed that the stationary states (i.e. MOs) are of the form
$$\psi = c_1\phi_1 + c_2\phi_2$$
where $\psi$ is an MO and $\phi_1, \phi_2$ 1s atomic orbitals on the two hydrogens, then one can use the variational principle to optimise the two coefficients $c_1, c_2$. The idea is that one wants to minimise the energy of the wavefunction $\psi$:
$$\mathcal{E} = \frac{\langle \psi | H | \psi \rangle}{\langle \psi | \psi \rangle}$$
with respect to both $c_1$ and $c_2$. To this end, both the following partial derivatives must be equal to zero:
$$\left(\frac{\partial \mathcal{E}}{\partial c_1}\right)_{\!c_2} = \left(\frac{\partial \mathcal{E}}{\partial c_2}\right)_{\!c_1} = 0$$
The maths will be omitted and if one is interested, one may find it in a variety of textbooks. For example it is covered in Atkins & Friedman's Molecular Quantum Mechanics 5th ed. (search for Rayleigh–Ritz method), or in Coulson's Valence 3rd ed. It is also described thoroughly in this handout from QMUL. The key result is the secular equations, which are a set of generalised eigenvector equations:
$$\mathbf{Hc} - \mathcal{E}\mathbf{Sc} = \mathbf{0}$$
where the $2 \times 1$ vector $\mathbf{c}$ simply contains the coefficients $c_1$ and $c_2$, the elements of the Hamiltonian matrix are $\mathbf{H}_{ab} = \langle a | H | b \rangle$ and the elements of the overlap matrix are $\mathbf{S}_{ab} = \langle a | b \rangle$. The basis set used is $\{|\phi_1\rangle, |\phi_2\rangle\}$. For a non-trivial solution (i.e. $\mathbf{c} \neq \mathbf{0}$) we require the secular determinant to vanish:
$$|\mathbf{H} - \mathcal{E}\mathbf{S}| = 0$$
Now we can use symbols to fill in these matrix elements:
The exact form of the Hamiltonian, as well as more explicit expressions for these quantities, are given in Atkins MQM. For now they need not concern us. Solving the determinant we obtain two values of $\mathcal{E}$ and two associated vectors $\mathbf{c}$, which can simply be written as the MOs $\psi = c_1\phi_1 + c_2\phi_2$. The MOs (i.e. approximate stationary states) are now quoted without proof, with their eigenvalues (i.e. their approximate energies):
$$\begin{align} \psi_1 &= \frac{1}{\sqrt{2 + 2S}}(\phi_1 + \phi_2) & & & \left\langle \psi_1 \middle| H \middle| \psi_1 \right\rangle &= \mathcal{E}_1 = \frac{\alpha + \beta}{1 + S_{12}} \tag{1} \\ \psi_2 &= \frac{1}{\sqrt{2 - 2S}}(\phi_1 - \phi_2) & & & \left\langle \psi_2 \middle| H \middle| \psi_2 \right\rangle &= \mathcal{E}_2 = \frac{\alpha - \beta}{1 - S_{12}} \tag{2} \end{align}$$
$\mathcal{E}_1$ is the ground state energy (i.e. $\mathcal{E}_1 < \mathcal{E}_2$) since $\alpha$ and $\beta$ are both negative quantities (again MQM has more info).
Now we have actually reached the part where the original question can be answered. What is the stabilisation of the bonding MO? The energy of the constituent AOs is $H_{11} = H_{22} = \alpha$. So, the stabilisation must be:
$$\begin{align} E_{\mathrm{stab}} &= \alpha - \mathcal{E}_1 \\ &= \alpha - \frac{\alpha + \beta}{1 + S_{12}} \\ &= \frac{\alpha S_{12} - \beta}{1 + S_{12}} \end{align}$$
Similarly, the destabilisation of the antibonding MO is:
$$\begin{align} E_{\mathrm{destab}} &= \mathcal{E}_2 - \alpha \\ &= \frac{\alpha S_{12} - \beta}{1 - S_{12}} \end{align}$$
The quantity $S_{12}$ is positive (again it should be in MQM), so $E_{\mathrm{destab}} > E_{\mathrm{stab}}$.
The "incorrect" piece of information comes from a further simplifying assumption - namely, that $S_{12} \ll 1$. That is not the greatest assumption, since the values of overlap integrals can sometimes be quite non-negligible (for a $\ce{H-H}$ bond length of $74~\mathrm{pm}$, $S_{12} = 0.753$!), but it is a useful way to predict approximate energies and the general forms of the MOs.
If $S_{12} \ll 1$, then $1 \pm S_{12} \approx 1$ and equations $(1)$ and $(2)$ become:
$$\begin{align} \psi_1 &= \frac{1}{\sqrt{2}}(\phi_1 + \phi_2) & & & \left\langle \psi_1 \middle| H \middle| \psi_1 \right\rangle &= \mathcal{E}_1 = \alpha + \beta \tag{3} \\ \psi_2 &= \frac{1}{\sqrt{2}}(\phi_1 - \phi_2) & & & \left\langle \psi_2 \middle| H \middle| \psi_2 \right\rangle &= \mathcal{E}_2 = \alpha - \beta \tag{4} \end{align}$$
Now it turns out that, with the additional neglect of overlap approximation, $E_{\mathrm{stab}} = E_{\mathrm{destab}} = -\beta$. That's where the "destabilisation of antibonding MO = stabilisation of bonding MO" idea comes from.
The arguments above can be extended to the case where you have constituent AOs which are different. The main difference is that $H_{11} \neq H_{22}$, so you have to write $H_{11} = \alpha_1$ and $H_{22} = \alpha_2$. The result you obtain will be rather similar: the destabilisation of the antibonding MO is greater than the stabilisation of the bonding MO, but in the limit of $S_{12} \to 0$, they become equal.
