How does the proof of operator commutativity work with non-continuous operators? In some books, a proof that if two self-adjoint operators $A$ and $B$ share a common eigenbasis $\{\phi_n\}$, then they commute is given as follows :

For any $\phi_n$,
$$AB\  \phi_n = a_n\  b_n\ \phi_n  = BA\  \phi_n,$$
and so $AB - BA = 0$.
Then, for any vector $\psi = \sum\limits_{n=1}^\infty c_n\  \phi_n$,
$$(AB-BA)\ \psi =  (AB-BA)\sum\limits_{n=1}^\infty c_n\  \phi_n= \sum\limits_{n=1}^\infty \ (AB-BA)c_n \phi_n = 0$$

But how is the third step on the to last line valid, where the operator $C = (AB-BA)$ is moved across the summation sign? $C (\lim f_n) = \lim (C f_n)$ is only valid if $C$ is continuous. But here, while $C$ is linear and self-adjoint, it may not be continuous!
 A: It is possible to say something more precise than Martin's answer (that is correct however).
  The key-point is that self-adjoint operators are  closed operators.
  An operator $A: D(A) \to H$, with $D(A) \subset H$ a linear subspace of the Hilbert space $H$ is said to be closed if, for every sequence of vectors $f_n\in D(A)$ such that 
(1) $f_n \to f \in H$ as $n\to +\infty$  
and 
(2) $Af_n \to g \in H$ as $ n\to +\infty$ 
then $f\in D(A)$ and $Af=g$.
In other words if both $D(A) \supset \{f_n\}$ and $
\{Af_n\}$ converge, then we can swap the symbol of limit and $A$:
$$A \lim_{n \to +\infty} f_n =  \lim_{n \to +\infty} Af_n\:.$$
There are several other equivalent definitions of closed operators, but this is the more elementary one and it is self-explanatory.
It is clear that continuous operators are closed but the converse is generally false. 
Within the most elementary formulation of QM, observables are represented by self-adjoint (generally unbounded and non-globally defined) operators $A$ on a suitable Hilbert space. So $A=A^*$. From the definition of adjoint operator, it turns out to be obvious that $A^*$ is closed. Thus $A$ is closed as well.
This weaker version of continuity permits to justify several naive manipulations as those you mentioned.
Let us now suppose that A is self-adjoint (in that case $D(A)$ is dense) admitting a complete orthonormal system of eigenvectors $\phi_n$, where $A\phi_n = \lambda_n \phi_n$ (I am supposing that the space is separable just for the sake of simplicity, as what I am writing holds even dropping that hypothesis). 
Obviously $\phi_n \in D(A)$ and, it is possible to prove (for instance using the spectral theorem mentioned by Martin), that $\psi = \sum_n c_n \phi_n \in D(A)$ iff $\sum_n |c_n \lambda_n|^2 < +\infty$.
The latter inequality, in view of the very definition of Hilbertian basis, is equivalent to the fact that  $\sum_n c_n \lambda_n \phi_n$ converges to some vector $\phi\in H$.
Since $A$ is closed, one easily sees that:
$$ A\psi = A \left( \sum_n c_n \phi_n \right) = \sum_n c_n A\phi_n =
\sum_n \lambda_n c_n \phi_n = \phi\:.$$
We have proved that:
THEOREM. Let $H$ be a (complex separable) Hilbert space. Suppose that $A$ is a self-adjoint operator in $H$, generally unbounded and defined on a (dense) linear subspace $D(A)$ of $H$, admitting a complete orthonormal system of eigenvectors $\phi_n$, where $A\phi_n = \lambda_n \phi_n$.
If $\psi \in D(A)$ then:
$$A \psi = \sum_n \lambda_n \langle \phi_n|\psi \rangle \phi_n\:.$$
Now, applying twice the theorem, we have an immediate corollary remembering that, for a pair of operators $C,D$ on $H$ with domains $D(C)$ and $D(D)$ respectively:
$$D(CD):= \{ f \in D(D) \:|\: D(f) \in D(C)\}\:.$$
COROLLARY. Referring to the theorem above, suppose that $B: D(B) \to H$ is another self-adjoint operator defined on the dense linear space $D(B) \subset H$ such that $B\phi_n = \mu_n \phi_n$ for every $n$.
If $\psi \in D(AB) \cap D(BA)$ then:
$$(AB-BA) \psi =0\:.$$
A: Since I'm not an expert on spectral theory, this will only be a partial answer, however, I believe that this question, is mathematically much more involved than you think.
First of all, let's review the finite dimensional case: We have two Hermitian matrices $A,B\in\mathcal{M}_d$ and they commute if and only if their spectral projections commute, i.e. they have a common eigenbasis. So much, so good.
Now we move to a separable, complex Hilbert space $\mathcal{H}$. The first thing to note is that for general bounded operators, the term "eigenvectorbasis" is not well-defined anymore. A bounded operator $A\in\mathcal{B}(\mathcal{H})$ might have continuous spectrum. So let's do something simpler first:
Definition: An operator $C\in\mathcal{B}(\mathcal{H})$ (bounded operator) is said to be compact, if it can be approximated by finite rank operator (i.e. is the "limit" of finite dimensional operators). Furthermore, a compact operator is called self-adjoint if it is Hermitian (i.e. $\langle \psi, H \phi \rangle = \langle H\psi, \phi \rangle$ for all vectors in the Hilbert space).
There are more natural equivalent characterizations, but I don't care about them here. Now we have the spectral theorem:
Theorem: Given a self-adjoint, compact operator $A\in\mathcal{B}{(\mathcal{H})}$, there exists a sequence of real numbers $\{\lambda_i\}_{i\in\mathbb{N}}$ accumulating at zero and a sequence of spectral projections $\{P_i\}_{i\in\mathbb{N}}$ such that 
$$ A=\sum_{i=1}^\infty \lambda_i P_i$$
Now, given this theorem, we can always and nicely use the proof above. We just need to note that any linear, bounded operator is automatically continuous (in fact, a linear operator is bounded iff it is continuous). 
Now let's get more complicated. For arbitrary bounded operators (again: they are always continuous!), we have a different spectral theorem. In order to formulate this, we need the notion of a projection valued measure, which will give us the continuous version of the spectral theorem, incorporating those parts of the spectrum that are not eigenvalues and therefore do not have eigenvectors (an example of such continuous spectrum can be seen for the Laplacian: Its spectrum is the positive part of the real line and the "eigenfunctions" would be non-square-integrable functions - $e^{ikx}$). So we define:
Definition: A projection valued measure on $\mathbb{R}$ is a map $\mu: \mathfrak{B}\to \mathcal{B}(\mathcal{H})$, where $\mathfrak{B}$ is the Borel sigma-algebra over the reals, which fulfills the following:


*

*$\mu(\mathbb{R})=1$ 

*for every $\xi,\zeta\in\mathcal{H}$, $\langle \xi, \mu(.)\zeta\rangle$ is a complex-valued measure.


Then we can formulate the theorem:
Theorem: Given a self-adjoint bounded, linear operator $A\in\mathcal{B}(\mathcal{H})$, there exists a projection valued measure $P_{\lambda}$ on $\mathfrak{B}$ such that 
$$ A=\int_{\mathbb{R}} \lambda dP_{\lambda} $$
This is sort of the continuous analogue of writing the operator as a sum of eigenvalues and projections. In the case of compact operators, our first theorem just tells us that the measure $dP$ is discrete. Now, a decomposition of a Hilbert space function into "eigenvectors" becomes less obvious - I don't see a really direct analogue to your proof above. However, we can do something similar. Note that the quantity $[A,B]$ is still well defined for every operator $A,B\in\mathcal{B}(\mathcal{H})$, so it makes sense to ask whether or not $[A,B]=0$. It seems that you can prove the following (the proof follows "easily" from the spectral theorem):
Proposition: Given two self-adjoint, bounded linear operators $A,B$, their spectral measures commute iff the operators commute.
Why do we want commuting measures? Well, in this case, it seems, one should be able to find a joint projection valued measure for the two observables! At this point, I'm a bit unsure of how to formulate this. I suppose, you should define the joint measure to be
$\mu: \mathfrak{B}\times \mathfrak{B}\to\mathcal{B}(\mathcal{H})$ is a projection valued measure such that $\int_{\mathbb{R}} d\mu(\lambda)=P_A$, where $P_A$ is the measure belonging to $A$ and similarlry, if I integrate over the other variable, I get the other measure. I also guess that the measure in question here should be just the product measure.
This should be as far as we can go with the analogy. Now, you were specifically asking for noncontinuous operators. In this case, we need to consider unbounded operators - these occur naturally as Hamiltonians in quantum mechanics (however, you can practically avoid them by always considering time evolution). First a definition again:
Definition: An linear operator $A:\mathcal{D}\to \mathcal{H}$ with $\mathcal{D}\subset \mathcal{H}$ is said to be unbounded, if for any $M>0$ there exists a $\psi\in\mathcal{D}$ such that $\|A\psi\|>M$. An unbounded operator is called Hermitian, if for all $\psi,\phi\in\mathcal{D}$ we have $\langle \psi , A\phi \rangle = \langle A\psi, \phi \rangle$. An operator is called self-adjoint, if in addition, the domain of its adjoint operator is also $\mathcal{D}$.
That is sort of clear, but I wanted to write this definition down to make clear the biggest difference between unbounded and bounded operators: While for bounded operators, we can always assume that the domain is the whole Hilbert space, this is no longer true for unbounded operators. In fact, unbounded operators are only defined on a dense subset $\mathcal{D}$ of the whole Hilbert space. At this point, operator multiplication becomes problematic - a quantity like $[A,B]$ is not a priori well-defined. You can of course easily make it a well-defined expression by just saying that the domain of $[A,B]$ is all vectors where it makes sense. However, in many cases, this will mean that the domain of $[A,B]$ is just $\{0\}$, which makes the expression $[A,B]=0$ a bit meaningless.
It turns out that you can use the same spectral theorem as above for unbounded operators (the projection valued measure for bounded operators will have compact support, while for unbounded operators, this is no longer true). Now we can have a different definition of commuting (this is now taken from Reed-Simon VIII.5):
Definition: Two self-adjoint operators $A, B$ commute iff all their projections in their associated projection-valued measures commute.
And then you can prove the following theorem (again Reed-Simon):
Theorem: $A,B$ two self-adjoint operators commute iff their one-parameter unitary groups commute (i.e. $e^{itA}e^{isB}=e^{isB}e^{itA}$ for all $s,t$). 
Given this, it might be reasonable to ask the question of whether there exists something like a "joint measure", but this might just as well be unreasonable. This is more than my knowledge and I guess you'd have to ask the mathematicians, if you wanted a more lucid and complete exposition of these subtle points...
As already stated, the last part is taken from
Michael Reed, Barry Simon: Methods of Modern Mathematical Physics I: Functional Analysis (chapter VIII).
The rest is from my head and might be faulty, comments and suggestions are very welcome. I can always provide more references and proofs, however most proofs (especially the spectral theorem) are quite complicated.
EDIT: In general, given a calculation of the type $A\sum_n a_n |f_n\rangle = \sum_n a_n A|f_n\rangle$, this can only be done for any state if $A$ is bounded. You will encounter this situation mostly in two contexts: either the operator is really bounded (then mostly it will even be compact and the $|f_n\rangle$ are the eigenstates) or the operator is unbounded and the $|f_n\rangle$ will be generalized eigenstates. Since these are not in the Hilbert space anyway, the whole expression is only formal and to see how this works mathematically, you have to work differently. 
Now there are situations different than this. One that I can think of is the Hamiltonian for the harmonic oscillator, which is unbounded, but has only pure point spectrum. In this case, for general coefficients, $H \sum_n a_n |f_n\rangle=\sum_n a_n H|f_n\rangle$ is mostly nonsense in the sense that the right hand side is infinite and the left hand side is not defined, but choosing the coefficients to be decaying fast enough, the right hand side makes actually sense. Then, you can easily show that it is a Cauchy sequence, hence by completeness converges to the left hand side. This however, has to be shown for any particular example of a sequence. In this particular case, you can still say that $H \sum_n a_n |f_n\rangle=\sum_n a_n H|f_n\rangle$ holds by just defining the relation to hold true if the right hand side is infinite and keeping in mind that this just means that the state is not in the support of $H$. However, I don't think this is general in any way.
