Modular forms could be defined for arbitrary subgroups of the modular group, and I have read that this is done in some papers, but every definition of a modular form I have seen has been with respect to congruence subgroups.

At a certain level, it's mostly a matter of (i) terminology and (ii) reading the right books. Technically the word "modular" in modular forms refers to the "modular group $SL_2(\mathbb{Z})$". In Miyake's book Modular Forms, he defines an automorphic form with respect to an arbitrary Fuchsian group $\Gamma$ (i.e., a discrete subgroup of $SL_2(\mathbb{R})$). Then he goes on to say (p. 114) that "Automorphic functions and forms for modular groups are called modular functions and modular forms respectively." Despite the title, plenty of the book deals with the general case, or with the special case of Fuchsian groups associated to quaternion algebras, which do not yield modular forms according to his definition. In Shimura's book Introduction to the Arithmetic Theory of Automorphic Functions he defines (pp. 2829) automorphic functions and forms with respect to an arbitrary Fuchsian group of the first kind (i.e., finite hyperbolic covolume). The phrase "modular forms" is sometimes used in his book, but doesn't appear to get a formal definition. These are, to my mind, the two most standard and authoritative references on "modular forms", and they both entertain the concept of a modular form with respect to a rather general Fuchsian group, whatever they want to call it. On the other hand, there are reasons for restricting to Fuchsian groups which are arithmetic (which is a technical term here) and of congruence type. A theorem of Margulis shows that arithmeticity is equivalent to having a sufficiently rich theory of Hecke operators, which is highly important in numbertheoretic applications. Similarly, being arithmetic of congruence type puts you in the realm of Shimura varieties, and gives you a rich theory of models of the Riemann surfaces defined over various abelian number fields. On the other hand, it is a famous consequence of Belyi's theorem that every algebraic curve over $\mathbb{Q}$ can be uniformized by a finiteindex subgroup of $SL_2(\mathbb{Z})$ (generally of noncongruence type). So if one is interested in the "special" arithmetic properties of modular curves, it makes sense to restrict to congruence type. Indeed, continuing work of John Voight and myself indicates that the congruence type condition is even more arithmetically significant than the arithmeticity [sic!]. We define congruence subgroups of nonarithmetic Fuchsian triangle groups and derive some of the arithmetic applications (using techniques from group theory and the arithmetic theory of branched coverings) that are parallel to those satisfied by the usual modular curves. See http://math.uga.edu/~pete/triangle091309.pdf Note that this work is not yet finished, to my consternation. (Mea culpa. Mea culpa.) 


The main point is that the basic definitions work fine but the link with arithmetic is much more "vague". Look at early papers of Tony Scholl. There are Galois representations attached to certain noncongruence forms, but they will in general only be of a subgroup of the absolute group of Q and they're typically not 2dimensional. So in summary, you can study them, sure, as Scholl and others did, but there are fewer applications. Here's an outline of what goes wrong. Modular curves defined by congruence subgroups have natural models over Q. A theorem of Eichler and Shimura says that the Tate module of the Jacobian of such a curve breaks up, as a Galois module, into "2dimensional bits". In the noncongruence case the curves aren't defined over Q and the galois representation on the Tate module doesn't break up into 2d pieces. 


Actually, modular curves for congruence groups have canonical models over number fields, not Q (there exist congruence subgroups other than $\Gamma _0(N)$!). They even have reasonably nice integral models. Moreover, the modular forms are sections of a sheaf defined on the canonical model, and the sheaf extends to the integral model. This has the following consequence: if a modular form has Fourier coefficients $a_n$ in a number field $K$ (so the form is a section of the sheaf over $K$), then the $a_{n}$ have bounded denominators, i.e., lie in $d^{1}\mathcal{O}_{K}$ for some $d$. Modular curves for arithmetic groups also have models over number fields, but the last statement definitely fails for forms that don't come from congruence groups. So something goes wrong with the beautiful picture we have for congruence modular curves, but I've never understood exactly what. However, this is another indication that the link to arithmetic is more tenuous in the noncongruence case, and helps explain why number theorists are mainly interested in the congruence case. [This was written as a comment on Buzzard's answer, but the site wouldn't let me post it (too long).] 


Another reason it is relatively uncommon to study automorphic cusp forms for noncongruence subgroups is that their existence is difficult if not impossible to establish, from the work of Phillips and Sarnak and others. EDIT: In the above I was referring to functions on the upper half plane invariant under the group $\Gamma$ which are eigenfunctions of the hyperbolic Laplacian, i.e., Maass forms. I should have been more explicit! For the case of $\Gamma = SL_2(\mathbb{Z})$ their existence was only shown by Selberg using the trace formula. The trace formula holds for noncongruence subgroups too but in general the spectral side has both a discrete part (corresponding to cusp forms) and a continuous part (corresponding to Eisenstein series). For $\Gamma_0(N)$ one can show directly that the continuous part contributes a lower order of magnitude than the whole spectral side, and therefore deduce that Maass forms exist in abundance. For noncongruence subgroups this argument breaks down and the work referred to above provides evidence that the discrete spectrum is very small. 

