Why is addition of observables in quantum mechanics commutative? I am no expert in the field. I hope the question is suitable for MO.
Background/Motivation
I once followed a quantum mechanics course aimed at mathematicians. Instead of the usual motivations coming from experiment at the turn of the 19th century, the following argument (more or less) was given to show that the QM formalism is in some sense unavoidable.
When one does physics, he is interested in measuring some quantity on a given state of the universe. The quantity (say the speed of a particle) is defined experimentally by the tool used to do the measure. We define such an instrument, with a given measure unit, an observable. So for every state and every observable we get a real number.
We can now define a sum and a product of observables. These are obtained by performing the two measures and then adding or multiplying their values. Similarly we can define scalar multiplication. These operations are then associative, but there is no reason why they should be commutative, since performing the first measure can (and indeed does) change the state of the universe. For some reason I cannot understand, anyway, addition is assumed commutative. I also see no reason why multiplication should distribute over addition. We can now also consider observables with complex values, by linearity.
At this point observables form an $\mathbb{R}$-algebra. We intoduce a norm it as follows. The norm of an observable is the sup of the absolute values of the quantities which can be measured. Every instrument will have a limited scale, so this is a real number. By definition this is a norm. Moreover it satisfies $\|A B \| \leq \|A\| \| B \|$. We can now formally take the completion of our algebra and obtain a Banach algebra.
Finally we define an involution * on our algebra by complex conjugation of observables. This yields a Banach * -algebra, and the third assumption which is mysterious to me is that the $C^*$ identity holds.
Finally we can use the Gelfand-Naimark theorem to represent the given algebra as an algebra of operators on a Hilbert space. If this turns out to be separable, it is isomorphic to $L^2(\mathbb{R}^3)$ and we recover the classical Schrodinger formalism.
The problems
In this approach I see three deductions which seems arbitrary: addition is commutative, multiplication is distributive and the $C^*$ identity holds. Is there any kind of hand-waving which can jusify these? In particular

Why is addition of observables commutative, while multiplication is not?

 A: This story cannot possibly be correct as written, and I agree with Fabian Besnard about why:

When you say "We can now define a sum and a product of observables. These are obtained by performing the two measures and then adding or multiplying their values."


This cannot possibly describe the usual sum A+B and product AB of operators.

Even simpler than what Fabien pointed out, this description if taken seriously would imply that the product of self-adjoint operators remains self-adjoint which is of course not at all true if they don't commute.
I actually don't believe a story of this form can be used to motivate the study of $C^{\ast}$-algebras, because it is not at all clear what the physical meaning of either addition or multiplication of observables is. It is not about adding or multiplying measurements because if two operators don't commute then their eigenvalues neither add nor multiply.
Let me propose an alternative (although I won't quite be able to finish it): instead of talking about measurements let's talk about Noether's theorem. Noether's theorem asserts, roughly, that to every $1$-parameter group of symmetries $\varphi_t$ of a classical mechanical system one can associate an observable $A$ such that the Poisson bracket $\{ A, - \}$ generates $\varphi$, and which is conserved, or equivalently, which Poisson-commutes with the Hamiltonian. Moreover, if the Hamiltonian of the system is time-independent then applying Noether's theorem to time translation symmetry produces the Hamiltonian $H$ itself, and we recover the fact that the Poisson bracket $\{ H, - \}$ generates time evolution. Similarly by applying space translation or rotation symmetry we recover conservation of momentum or angular momentum respectively.
Now we can ask: what other mathematical objects satisfy a form of Noether's theorem?
As motivation consider the following: the Skolem-Noether theorem implies that the simple matrix algebra $M_n(\mathbb{C})$ has the property that all of its automorphisms as a $\mathbb{C}$-algebra are inner. This means a one-parameter group $\varphi_t$ of automorphisms must be a one-parameter subgroup of $PGL_n(\mathbb{C})$, hence is generated by some element of the Lie algebra $\mathfrak{pgl}_n(\mathbb{C})$. By lifting this element to $\mathfrak{gl}_n(\mathbb{C})$ we get that there must exist an $A \in M_n(\mathbb{C}) \cong \mathfrak{gl}_n(\mathbb{C})$ such that
$$\varphi_t(X) = \exp(At) X \exp(-At).$$
Now, any one-parameter group of automorphisms of a finite-dimensional algebra must be generated by a derivation, and differentiation shows that the derivation generating $\varphi_t$ is the commutator $[A, -]$. In other words, Noether's theorem holds for $M_n(\mathbb{C})$ with the Poisson bracket replaced by the commutator! We can even get the Hamiltonian into the game by positing that it is just some other element $H \in M_n(\mathbb{C})$ such that the commutator $[H, -]$, as above, generates time evolution, meaning that time evolution is a one-parameter group of automorphisms of $M_n(\mathbb{C})$, and then we find that if some other one-parameter group of automorphisms preserves $H$ then it corresponds to some $A$ which commutes with $H$: $[A, H] = 0$.
This is, to my mind, an extremely suggestive result, for a state of mind in which we know the Hamiltonian formulation of classical mechanics but know nothing whatsoever about quantum mechanics. It tells us that $M_n(\mathbb{C})$ and classical mechanics, while outwardly looking very different, have the following features in common:

*

*A Lie bracket, either the commutator or the Poisson bracket, acting on some vector space (of "observables"?).

*A notion of one-parameter groups of automorphisms of the space of observables, which can be generated by the Lie bracket in the sense that the brackets $[A, -]$ or $\{ A, - \}$ produce derivations which can be exponentiated to automorphisms. In particular time evolution is itself such a one-parameter group of automorphisms and the corresponding generator is, in classical mechanics, the Hamiltonian $H$.

*A Noether's theorem, in the following two-part form: every one-parameter group of automorphisms comes from the Lie bracket, and two such groups commute with each other iff the corresponding bracket of their generators is zero.

To my mind, what this line of reasoning suggests is that the really meaningful operation on observables is not an associative multiplication but a Lie bracket, not because it has a direct interpretation in terms of measurements but because a Lie bracket is the sort of thing you need to state a Noether's theorem, and in particular to convert a specific observable $H$ called the Hamiltonian into a one-parameter group of automorphisms describing time evolution.
This suggests a potential route to an answer to the title question about why addition of observables in quantum mechanics must be commutative: it's because addition in a Lie algebra must be commutative (I see that Gianmarco Bramanti has also suggested this). I like this idea because in the Lie algebra case it's also unclear what the "physical" meaning of addition is: it doesn't correspond in a straightforward way to an operation on one-parameter groups if the corresponding generators don't bracket to zero, which is exactly the issue we see with observables.
So, why is the Skolem-Noether theorem true anyway? The proof is short and elegant: $M_n(\mathbb{C})$ has a unique simple module, namely $\mathbb{C}^n$. Therefore if $\varphi : M_n(\mathbb{C}) \to M_n(\mathbb{C})$ is any automorphism, then pulling $\mathbb{C}^n$ back along $\varphi$ produces a module which must be isomorphic to $\mathbb{C}^n$ again. If we require $\varphi$ to be $\mathbb{C}$-linear then, since $M_n(\mathbb{C})$ is the full algebra of $\mathbb{C}$-linear endomorphisms of $\mathbb{C}^n$, this isomorphism must itself be given by an invertible element of $M_n(\mathbb{C})$, and then we can check that this means $\varphi$ must be conjugation by this element.
This suggests that we try looking for infinite-dimensional algebras like $M_n(\mathbb{C})$ which also have unique modules in some sense, which are also the full algebra of endomorphisms of those modules in some sense, and which have enough analytic structure to allow us to talk about exponentiation of derivations. We need algebras of endomorphisms specifically to run the second half of the argument above.
Now (this is the unfinished part, I don't have a clear sense of how to tell the story from here) we can bring in algebras of operators on Hilbert spaces, the Stone-von Neumann theorem about representations of the canonical commutation relations $[X, P] = i \hbar$, Stone's theorem on one-parameter unitary groups, and the Skolem-Noether theorem for the algebra of bounded operators on a Hilbert space, which together are also quite suggestive, although the leap to Hilbert spaces and the self-adjoint / unitary constraints haven't been motivated (for that matter we haven't at all motivated the decision to work over $\mathbb{C}$ as opposed to $\mathbb{R}$). This still does not get us to $C^{\ast}$-algebras but IMO this is a feature and not a bug: $C^{\ast}$-algebras can only talk directly about bounded operators and the position and momentum operators aren't bounded! Also there's no reason an arbitrary $C^{\ast}$-algebra should satisfy Noether's theorem; this line of reasoning specifically pulls us towards the special ones that do, which IMO is also a feature and not a bug.
A: We may think of a state $\omega$ as a functional on the algebra of observables $\mathcal O$ which is interpreted as giving the expected value of each observable. With this in mind, it is natural to require $\omega$ to be linear (as well as two other usual properties, positivity and normalization).
Thus given two observables $A, B \in \mathcal O$, if we are going to have a sum $A + B$ it should be true that $\omega(A + B) = \omega(A) + \omega(B)$ for any state $\omega$. Since this gives the values of $A + B$ on every state, it suffices to define it. Since $B + A$ has the same values on every state, $B + A$ is the same observable.
On the other hand, there is no natural way to say what $\omega(AB)$ should be, since states (like expectation values) need not be multiplicative.
So: commutativity of observables reduces to commutativity of $\mathbb C$ since expectations are linear, but nothing analogous applies to multiplication.
This is based on what I've read in F. Strocchi, An Introduction to the Mathematical Structure of Quantum Mechanics.
Note that you can eventually interpret states as arising out of probability distributions, leading to Theo's comments.
Personally I am still a bit hazy on why we postulate a multiplication on observables at all, when (unlike the classical case) there is not a clear physical interpretation of what such an operation should mean. However, given the full structure of a $C^*$-algebra, one can show that the uncertainty principle (or the existence of complementary observables) requires noncommutative multiplication, et voila, you have quantum mechanics.
A: This question has bothered me for a long time! Although I don't have an answer, I'd like to mention an approach that looks promising at first, but turns out not to work.
First, recall that in quantum mechanics, you can think of a "state" as a way of preparing a physical system. Theo Johnson-Freyd pointed out in a comment that if you have two states $\rho$ and $\sigma$, you can construct a state that intuitively deserves to be called $\tfrac{1}{2}(\rho + \sigma)$:

Flip a fair coin. If the coin comes up heads, prepare the system in state $\rho$. If the coin comes up tails, prepare the system in state $\sigma$.

This state deserves the name $\tfrac{1}{2}(\rho + \sigma)$ because if $\rho[X]$ is the expectation value of the observable $X$ for a system prepared in state $\rho$, and $\sigma[X]$ is the expectation value of $X$ for a system prepared in state $\sigma$, the expectation value of $X$ for a system prepared in state $\tfrac{1}{2}(\rho + \sigma)$ should be $\tfrac{1}{2}(\rho[X] + \sigma[X])$, by the laws of classical probability.

Now, what happens if we use the same trick to define the sum of two observables? Given two observables $X$ and $Y$, let's define $X + Y$ to be the observable:

Flip a fair coin. If the coin comes up heads, measure $X$ and double the result. If the coin comes up tails, measure $Y$ and double the result.

The laws of classical probability tell us that if $\rho[X]$ and $\rho[Y]$ are the expectation values of $X$ and $Y$ for a system prepared in state $\rho$, the expectation value of $X + Y$ for a system prepared in state $\rho$ should be $\rho[X] + \rho[Y]$, just as you would hope.

Here's where things go pear-shaped. Given an observable $Z$, it makes sense to define $Z^2$ to be the observable:

Measure $Z$ and square the result.

So what's the expectation value of $(X + Y)^2$ for a system prepared in state $\rho$? The laws of classical probability tell us that it's $\rho[X^2] + \rho[Y^2]$. In the formalism of quantum mechanics, however, $X$ and $Y$ are operators and $\rho$ is a linear functional on the operator space, so

$\rho[(X + Y)^2] = \rho[X^2] + \rho[Y^2] + \rho[XY + YX]$.

If $\rho[XY + YX]$ is nonzero, this formula disagrees with the expectation value for $(X + Y)^2$ that follows from our definitions of $X + Y$ and $Z^2$, according to the laws of classical probability!
In practice, it's not hard to find observables $X$ and $Y$ for which $\rho[XY + YX]$ can be nonzero. For example, let $X$ and $Y$ be the x-spin and z-spin of a spin-1 particle, represented by the operators
$X = \frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0&1&0\\\\1&0&1\\\\0&1&0\end{array}\right],\qquad Y = \left[\begin{array}{ccc}1&0&0\\\\0&0&0\\\\0&0&-1\end{array}\right].$
A: If I correctly understand, you are misinterpreting the meaning of the product and sum of observables.
When you say "We can now define a sum and a product of observables. These are obtained by performing the two measures and then adding or multiplying their values."
This cannot possibly describe the usual sum A+B and product AB of operators. For the product, it is not even hermitian unless A and B commute. Agreed, A+B is hermitian, but the spectrum of A+B does not contain the result of the sum of a measurement of A followed by a measurement of B (in either way), again unless A and B commute. For a counter-example take $A=\pmatrix{1& 0\cr 0&-1}$ and $B=\pmatrix{0&1\cr  1&0}$.
I hope I correctly understood your question.
A: The introduction you outlined is basically reminiscent of the old quantum mechanics, anyway in the approach you depicted it is the culmination and not the premise of the construction, and the comment was surely intended to be explanatory of the far origin of these choice. I now try to resume the history.
There are basically two approaches to mathematical quantum mechanics. The first one very complex and stratified in its development, but simple in the premise was discussed by John Von Neumann in a lot of papers after the "Foundation of quantum mechanics", the second one is basically conceveid to be an extension for the first, and is this second approach you are referring to: the GNS approach. 
Anyway both of them are surely derived after an abstraction process very far beginning on the methods of classical mechanics, joint to the newest evidence from atomic and particle physics of the first quantum mechanics.  
Just as in classical mechanics we define functions of observables dynamical quantity so the founders of quantum mechanics conceived it is possible in quantum mechanics, anyway we need to clarify in which sense this is possible and explanation isn't fully depleted from the naive extension of classical theory of the measure, based on real numbers, but it need of a clear axiomatic and this was furnished from John Von Neumann (and in some way from Heisenberg, Dirac, and Schroedinger before him formulated this axiomatic)
Anyway, just as in classical mechanics there is a notion of repeatibility and regularity, so there is in quantum mechanics. The true difference is in the outcome of the measures, deterministic in classical, probabilistic in quantum mechanics. So that measure processes are conceived deterministic in a statistical sense, and, for example, the component energies of the isotropic harmonic oscillator sums exactly in mean value, but the variance is zero if the considered states are eigenstates. Old quantum mechanics can be founded on few  axioms about the measures and led Von Neumann, in a natural way to linear operators acting, like a non commutative algebra, on Hilbert spaces. 
In order to grant correspondence principle we, following the founders of quantum mechanics, need to hypothesize the existence of intrinsically deterministically evolving  observable, and just the measure process make the difference, because these dynamical "quantities" with respect to the measures doesn't appear as real numbers, this point was the first time realized some time after the Copenaghen interpretation was developed.
So they are assumed, after Heisenberg (speaking of non commutative numbers)  and Jordan (speaking of matrices), and Schroedinger (speaking of operators acting on functional space of probability) all these three point of view were showed to be in a certain strict framework to be equivalent, from Dirac assuming they are algebraic elements obeying to canonical commutation relation generalizing the Poisson algebra. 
In brief the Dirac point can be summarized in assuming an Hilbert space structure for the states, and in developing step by step a theory of observables compatible with the Copenaghen interpretation spirit and with the correspondence principle. 
Anyway Von Neumann felt the need to obtain an axiomatic foundation  based on more general operators algebras, and an axiomatic of measure, unifying from scratch the theoretical 
framework, in  fact obtaining a more general theory with respect to the Heisenberg and Dirac theoretical "prejudices". The Von Neumann point was in fact based on the general representation theory in the geometrical framework of Banach operator algebras of operators in Hilbert space, and in particular on the CCR irreducible representation theory, but from this point the research of Von Neumann continued in search of an intrinsic point of view based on the geometry of observable.
After time and time was in fact recognized that part of quantum theory of measure is nothing else then a generalized probabilistic theory in a Banach algebra and the general setting of Gelfand Najmark Segal construction rebuild intrinsically the Hilbert spaces. Anyway the field extension of this setting is very problematic and a hierarchy of Hilbert spaces appears. Anyway in this way a circle is closed and a new loop is opened: in the GNS approach to quantum mechanics we postulate that operators are living in an abstract algebra, obeying familiar rules for an algebra with an involution (the * operation). Via Gelfand theorem the commutative case led to the algebra of complex valued continuous functions in an Hausdorf space, the spectrum of the algebra (which will led the ordinary numerical set of coordinates of classical mechanics), and more in general to a spectral theory, culminating in the GNS construction, which associate to a given linear form an Hilbert space and a representation for the algebra. 
Anyway the true achievement of this approach is the net of algebras, that is very more general with respect to the Hilbert space interpretation of quantum mechanics,this achievement is useful in relativistic field theory and leads to very far reaching results firstly partially discovered  by Von Neumann in some papers, and after then developed from Araky, Haag, Kastler.  In this full setting is now possible to address in more precise terms the question of the cluster decomposition principle implicit in the deterministic evolutionary scheme, and the question of repeatability principle of classical and quantum mechanics, and to understand quantitatively something about the limitation, arising from the change of the state of the universe, to this principle, which can be espressed, for example, in term of a change of representation, becaused from the change of the linear form representing the thermokinetic state of "universe", without any change in the postulates of quantum field theory and the derived quantum mechanics. This is perhaps  the perspective of the search about KMS theorem.
I'm not very satisfied from this resume, anyway I think you can correct and integrate it, and I hope to read and write something else more precise and delimited.   
A: Since you raised your question I'm uncomfortable about some question, and I re-read Von Neumann, in order to becalm myself, anyway in Von Neumann the problem isn't solved in a deductive way neither is justified at all, only is asserted the additivity, as customary, between commutating operators (so there the correspondence principle is granted) and then the additivity is extnded to non commutating operators. 
Anyway reflecting on the practical use of non commutating linear combination of operators I realize that the arguments of linear combination are generally elements of some Lie Algebra, and their powers, and I remember a point in the first book of Landau about classical mechanics that I like to quote:
Conservation laws.
"Not all integrals of motion have an equally relevant role in mechanics. Among these there are some whose invariance over time has an origin very deep, related to fundamental properties of space and time and that is their homogeneity and isotropy. These quantity, these conservative, have an important general property, they are additive, that is, their value for a system composed of several elements, whose interaction can be neglected, is equal to the sum of the values for each of the elements."
It seems just like Landau is mixign two unrelated points: the isotropy and the commutativity. In fact this isn't the additivity we are thinking to. Anyway there is an important point: and this is the role of simmetry and in mathematic the role of the Stone-Neumann Theorem, and the role of parallel transport, and in mathematic the role of gauging the space and time, so we can perahps  re-connect the two point arised by Landau to the additivity of observables. 
Generalizing, perhaps we need to relate an observable to an infinitesimal of a continous symmetry group: an elementar generator of lie algebra in order to justify additivity. Just handwaving.
What do you think about this?
A: Your description of the structure of the algebra of observables isn't quite how I'm used to it being.  Indeed, I believe that in the best algebraic descriptions of quantum mechanics, addition is a formal operation, rather than a physical operation as you've described.  The best reference I know for this point of view is L.D. Faddeev and O.A. Yakubovskii, 2009, Lectures on Quantum Mechanics for Mathematics Students.  I don't have my copy handy right now, so I will describe my memory of how they set up the algebra of observables.
The first thing to point out is that in the real world, there is no such thing as pure states.  This has nothing to do with quantum mechanics, and everything to do with an experimenter's inability to perfectly measure the initial set-up.  For your notion of "state" to make sense physically, it must be something like "repeatable initial set-up for an experiment".  Once this is your notion of state, you are perfectly able to run your experiment 1000 times, make your measurements (each individual run may give a different answer, but you can look at the distribution), and process them as you want.
So really an observable assigns a probability distribution on $\mathbb R$ to each state.  We now demand the following axiom: the (good) functions $\mathbb R \to \mathbb R$ act on the set of observables by composition.  So if $X: \{\text{states}\} \to \{\text{probability distributions}\}$ is an observable, so is $X^2$: the probability that the observable $X^2$ assigns to an interval $[a,b]$ is the same as the probability that $X$ assigns to the interval $[\sqrt a,\sqrt b]$.  In particular, suppose you compose your observable $X$ with a step function $\Theta(x - \xi)$, where $\xi \in \mathbb R$.  Then the observable $\Theta(X-\xi)$ measures the whether the value of $X$ is more than $\xi$.  Then you can check that the full distribution $X$ is recoverable from the knowledge of all the $\Theta(X-\xi)$.  In particular, it's recoverable from the expected values of $\Theta(X-\xi)$ on each state.  So to set up the algebra of observables, it's enough to know only the expectation values for observable at each state.
Now you should realize the following.  The previous paragraph makes sense even for classical mechanics, and in fact is the correct formalism (as there are no pure states).  But in quantum mechanics, it's worse than that.  A definite state is one that gives a delta distribution for each observable.  Classically, we believe that a sufficiently good experimenter can approximate definite states to whatever desired accuracy.  But there is very good evidence that this fails in the quantum world: no matter what tools you use, there are absolute bounds preventing states from approximating definite states.  So the language of distributions and expectations is absolutely necessary to formalize quantum mechanics, whereas in classical mechanics you could say that there are idealized definite states, observables are functions on definite states, and states are probability distributions on the space of definite states.
Finally, the question is how to assign algebraical operations to the collection of observables.  And here I admit that I don't have a great answer.  One possibility is simply to convolve probability distributions: this gives a commutative addition, for example.  Then you could define a commutative associative multiplication by taking logs and adding and exponentiation, but my memory is that this does not distribute over addition in general. F&Y define a commutative nonassociative multiplication by $(X,Y) = \frac12\bigl((X+Y)^2 - X^2 - Y^2\bigr)$.  Oh, right.  The problem is the following: do you add, multiply, etc. the distribitions, or the expectation values?  For addition, adding expectation values is the same as the usual convolution of distributions and then taking expectation.  But for multiplication it is not.  I don't remember what F&Y do, but I think it might at the level of expectation values.
A: Trying to understand the introduction of the quantum formalism I have come across this question myself. I paste here the answer that I have given to myself and that convinces me the most. Please, if there are errors, let me know so I can continue learning:
We start by assuming a classical setup:

*

*we have pure and mixed states (which are probability distributions of pure states).

*the observables are real-valued functions, which can be seen like part of a c-star algebra $\mathcal{A}$.
We can reverse our point of view and think of all this like if states were functionals on $\mathcal{A}$, specifically they are normalized, positive, linear functionals. That is, a mixed state $\omega$ is a probability distribution and an observable $f$ is a random variable, so we consider the pairing $\omega(f)=E_{\omega}(f)$ (the expected value). The pure states are those functionals which are also multiplicative (degenerated distributions, with all the probability concentrated on a single point). We decide to stick to this point of view, since it is equivalent to the original one (this is proven by the Gelfand-Naimark theorem and Riesz-Markov representation theorem), and because at the end of the day observables are what we manage in real life. If you can't measure the difference between two systems, you have no right to treat them as different. By the way, in mathematics there are lots of examples where a space is recovered from its algebra of functions, or at least from its sheaves. This is the spirit of Algebraic Geometry.

Now, some experimental facts, like for example the discrete emission spectrums, the Stern-Gerlach experiment, and so on, led people like Heisenberg to conclude that there were observables $A,B$ such that the mean square deviations $\Delta_{\omega}(A),\Delta_{\omega}(B)$ cannot be simultaneously reduced. For example, it was checked that for the position $q$ and momentum $p$ of electrons in any state $\omega$ it is satisfied the relation
$$
\Delta_{\omega} (q)\Delta_{\omega} (p)\geq h/4\pi \equiv\hbar/2, \tag{1}
$$
named Heisenberg uncertainty relation (think that when we watch an electron we modify its momentum). But this doesn't fit the description above, since in theory we could refine the apparatus and obtain a pure state $\omega$ such that $\Delta_{\omega} (q)=\Delta_{\omega} (p)=0$.
On the other hand, without having anything to do with this, observe the following algebraic development. Let $A=A^*$ and $B=B^*$ be two self-adjoint elements of an arbitrary c-star algebra such that, without lost of generality, satisfy $\omega(A)=\omega(B)=0$, being $\omega$ a pure state (normalized, positive, linear functional). For $\lambda \in \mathbb R$ we have $(A-i \lambda B)(A+i \lambda B)$ is a positive observable and positivity of $\omega$ implies
$$
\omega\left(A^2\right)+|\lambda|^2 \omega\left(B^2\right)+i \lambda \omega([A, B]) \geq 0,
$$
where we define, right now, $[A,B]:=AB-BA$. But if we think of this expression as a polynomial in $\lambda$ with real coefficients then we have:
$$
4 \omega\left(A^2\right) \omega\left(B^2\right) \geq|\omega(i[A, B])|^2,
$$
and then
$$
\Delta_\omega(A) \Delta_\omega(B) \geq \frac{1}{2}|\omega([A, B])|. \tag{2}
$$
So a possible explanation for the Heisenberg relation (1) could be that the observables do not belong to a commutative c-star algebra but to a non commutative one!!!! We can introduce the canonical commutation relations in the algebra to explain that empirical fact.
What is happening here?
As mathematicians (or as human beings), we elaborate theories assuming some things (and deducing others). At the beginning (classical physics), we chose modelling "experimental data outputs", or observables, as a commutative algebra: the functions defined on the phase space. This was the assumption. Since observables are tightly related to measurement apparatuses we thought that the operations of the algebra $\mathcal A$ should be given by the corresponding real number operations with the measured values.
But new experimental facts revealed that reality would be better modelled if we drop the "commutative product" hypothesis. We have that the observables "xposition" and "yposition" can be simultaneously measured, in real life, with precision enough, so we can let them behave "classically" from the operational point of view. But since "xposition" and "xmomentum" cannot be simultaneously measured with precision (as pointed out by Heisenberg (1)) then we should define multiplication in $\mathcal A$ in such a way that they do not commute, and then relation (2) would explain (1).
In other words, the product of observables cannot be the pointwise product of the obtained values, but some other operation which only when restricted to the particular case of compatible observables is the pointwise product. Does this product have an interpretation in terms of the output of the measurement apparatuses? I am not aware of any.
The sum remains being commutative because we don't have the necessity of remove this hypothesis, and we can still stick to the idea of observables living inside a big c-star algebra. In case of compatible observables like "xposition" and "yposition", the sum should be what you naturally would think: a new observable whose meaning is self evident. All the observables compatible with, for example, "xposition" constitutes a subalgebra (let's denote it by $\mathcal{A}_{xposition}$) where we can perform sums and products and still remain inside.
Now consider "xposition" a "xmomentum", for example. The sum is well defined even as a new observable, since the sum of self adjoint elements is self adjoint. But, what is its meaning?
To my understanding, the sum and product of observables belonging to "different subalgebras" is not naturally defined. It depends on the "big algebra" $(\mathcal A,+,\cdot)$ we choose. But keep an eye, this choice has to reconcile what we know from experiments:

*

*The sum and product on $\mathcal A$ are such that when restricted to any subalgebra $\mathcal A_B=C_{\mathcal A} (\{B\})$, where $B\in \mathcal A$ and $C$ denotes the centralizer, the sum and product of observables belonging to $\mathcal A_B$ corresponds to the sum and product of the outputs when they are interpreted as apparatuses.

*Some relations for elements belonging to different subalgebras, like $q\cdot p-p\cdot q=i\hbar$, have to be introduced "by hand", in order to algebraically justify Heisenberg uncertainty relations.

So the answer to "Why is addition of observables in quantum mechanics commutative?" is:
Because so far there has been no need to eliminate that hypothesis. But we cannot even be sure that we have more than several commutative algebras $\mathcal A_B$ tied together with some commutation relations. That is, not only can we not be sure that the sum of two arbitrary elements is commutative, but we cannot even be sure that such a sum exists. But as I understand it, at the moment there is no empirical contradiction in assuming it.
The following is a link to the same text but with hyperlinks that maybe help to follow the reasoning.
A: If you know how to measure $<A>$ , $<B>$ and $<A+B>$, then you can test by experiments, whether $<A+B> = <A> + <B>$ .
A nice example is the Stern-Gerlach experiment that you have to perform for 3 different directions.
