You are right to be questioning the validity of this method, and I congratulate you for doing so. This is actually an extremely important skill, that differentiates the best students from the rest.
The variational method, like many methods from physical chemistry, is a method of approximation (a model) to what really happens. There exist a whole hierarchy of methods for computing the molecular orbitals of molecules (which are themselves models, being the stationary states of the Schrödinger equation), including at the top end post-Hartree-Fock theory and density functional theory. These methods provide quantitative information about molecular orbitals, but also require serious computer power (which is generally unavailable to undergraduate students). However, all of the basic physics can be explored and understood with simpler models that can be solved on a few sheets of paper, such as LCAO theory, hence why we teach them to undergraduates.
As regards to your question, we do indeed optimise a linear combination. We start off by assuming that the total molecular orbital wavefunction can be approximated using a linear set of atomic orbitals:
$$|\Psi\rangle = c_1|\phi_1\rangle + c_2|\phi_2\rangle + \cdots + c_n|\phi_n\rangle$$
We then need to find the coefficients. The variational principle (also known as the Rayleigh-Ritz method) states that the coefficients that give the best approximation to the wavefunction will minimize the energy, given by
$$\mathcal{E} = \frac{\langle\Psi|\hat{H}|\Psi\rangle}{\langle\Psi|\Psi\rangle}$$
Now, the second part of your question involves how to compute the two terms in this fraction. Without going into massive mathematical detail, the LCAO method can be recast into a matrix problem rather than an integral problem using the atomic orbitals as a form of basis vector. In the simplest case (Hückel theory) we assume that they are normalized, such that the denominator is always 1 or 0 (see later). The problem now is mostly how to determine the numerator.
In brief, each element of the Hamiltonian matrix is given by
$$H_{ij} = \langle\phi_i|\hat{H}|\phi_j\rangle$$
The variational principle applied to this matrix implies that the optimized energies are the eigenvalues of the Hamiltonian matrix, and the coefficients are given by the eigenvectors. Since an $n$ x $n$ matrix has precisely $n$ eigenvalues, this implies the "number of molecular orbitals" rule. Notice that the reason for the +/- combinations is by our own design (LCAO). We could have picked more complex trial functions by assumption, but this would make computation far more difficult.
It most certainly is true that in order for there to be a significant interaction, two orbitals must be close in energy. The detailed reasons are complex, but essentially it comes down to the size of the $\langle\phi_i|\hat{H}|\phi_j\rangle$. Orbitals that are far apart in energy have small values of this term hence interact weakly. This doesn't necessarily mean they can't form molecular orbitals, however, but these effects are negligible and are not important in understanding the chemistry, so are neglected in simple models (although are often included in some of the most complex modern methods).
The reason that only atomic orbitals of the same symmetry give molecular orbitals is the overlap integral $\langle\phi|\phi\rangle$. This is the total sum (integral) of the product of the AOs. This is zero for different symmetries. This must be non-zero to give an MO. Consider a $\mathrm{2p_z}$-orbital and a $\mathrm{1s}$-orbital for instance (for illustration purposes only!). The s-orbital is spherically symmetric, with all points the same sign. The $\mathrm{p_z}$-orbital has a "dumbbell" shape, with equal areas of different signs above and below the $xy$-plane. Therefore, the total sum of the product is zero. No overlap = no interaction.
