We study the combinatorial problem which consists, given a system of linear relations, of finding a maximum feasible subsystem, that is a solution satisfying as many relations as possible. The computational complexity of this general problem, named MAX FLS, is investigated for the four types of relations =, greater than or equal to, > and not equal. Various constrained versions of MAX FLS, where a subset of relations must be satisfied or where the variables take bounded discrete values, are also considered. We establish the complexity of solving these problems optimally and, whenever they are intractable, we determine their degree of approximability. MAX FLS with =, greater than or equal to or > relations is NP-hard even when restricted to homogeneous systems with bipolar coefficients, whereas it can be solved in polynomial time for not equal relations with real coefficients. The various NP-hard versions of MAX FLS belong to different approximability classes depending on the type of relations and the additional constraints. We show that the range of approximability stretches from APX-complete problems which can be approximated within a constant but not within every constant unless P = NP, to NPO PB-complete ones that are as hard to approximate as all NP optimization problems with polynomially bounded objective functions. While MAX FLS with equations and integer coefficients cannot be approximated within p(epsilon) for some epsilon > 0, where p is the number of relations, the same problem over GF(q) for a prime q can be approximated within q but not within q(epsilon) for some epsilon > 0. MAX FLS with strict or nonstrict inequalities can be approximated within 2 but not within every constant factor. Our results also provide strong bounds on the approximability of two variants of MAX FLS with greater than or equal to and > relations that arise when training perceptrons, which are the building blocks of artificial neural networks, and when designing linear classifiers.