Energy of a Hamiltonian that is the sum of commuting Hamiltonians

If your Hamiltonian is a sum of commuting Hamiltonians as you write then it is guaranteed that there exists a basis that simultaneously diagonalizes all Hamiltonians together, and therefore $E=sum_{i} E_i$ indeed. If any of the $H_i$ is non-degenerate, then diagonalizing it will generate the unique basis of the eigenstates. However, they will not be $psi = prod_i psi_i$. I assume that in this you mean some tensor-product? If the case is that the parts of the Hamiltonian commute because they operate in orthogonal subspaces of the Hilbert space, then this will be the case. Then you will have to diagonalize each Hamiltonian and indeed proceed as you wrote.

Note, that if the Hamiltonians are degenerate, a choice that diagonalize one of them might not be diagonal in the other, just block-diagonal.

Clearly any state of the form $psi = otimes_i psi_i$ will be an eigenstate of $H$ with eigenvalue $E = sum_i E_i$. But it's possible that there are eigenstates of $H$ that don't factorize, in case there are degeneracies. Consider for instance the case where $H$ is the sum of two harmonic oscillators, having energies $omega_1 = 1$ and $omega_2 = 2$. Then for any $alpha,beta in mathbb{C}$ the state $$| psi rangle = (alpha (a_1^dagger)^2 + beta a_2^dagger) | text{vac} rangle $$ is an eigenstate with energy 2 (since $2omega_1 = omega_2 = 2$). (I'm assuming that $H$ is normal-ordered such that the ground state has energy $E = 0$.) Unless $alpha = 0$ or $beta = 0$, the state $|psi rangle$ does not factorize.

Properly speaking one should require that the spectral measures of single parts of the total Hamiltonian commute. This is however the standard if the parts correspond to independent subsistems, and I interpret this way your question: each part is defined in a single Hilbert space and the total space is the tensor product of them. What finally happens is that the spectrum of the total Hamiltonian is the closure of the sum of the spectra. I do not think one can say further. Restricting to the point spectrum, it turns out that sum of proper eigenvalues are eigenvalues, but you can obtain the same result summing different strings of eigenvalues in general. So your claim regarding the form of eigenvectors may fail. (The naive claim should however refer to a basis if eigenvectors of the complete Hamiltonian: in general, one has to take linear combinations of your candidates in case of degenerated eigenvalue of some component Hamiltonian.)

Referring to the general case, where the operators do not live in separated Hilbert spaces, the only result relies on this fact: if $H_1, ldots, H_n$ are selfadjoint and pairwise commute then there is a selfadjoint operator $A$ such that $H_i=f_i(A)$ for some (Borel) functions $f_i$. Hence $$H=sum_{i=1}^n f_i(A) =: f(A).$$ Hence the spectrum of $H$ is the closure of the image of $R$ through $f$, that is the closure of the sum of the spectra of the $H_i$.