If he were alive today, Sylvester might have tweeted: The quaternion equation ax+xb=c is the simplest nontrivial linear equation in this noncommutative division ring. It provides a "canary in the coal mine" moment for newbies to noncommutative algebra. Yes, we can translate this equation into a 4D matrix problem, and utilize all of the tools of matrix algebra to solve it, but it doesn't seem very satisfying to not solve the equation *within* the noncommutative algebra of quaternions. The problem with Sylvester's equation is that since we can't get both the a and b coefficients on the same side of the unknown x, we can't utilize the distributive law to put the unknown x into a single term. But by using an extremely clever feline trick, we can group the terms in such a way that we can commute *scalar* numbers such as N(a)=aa', N(b)=bb', S(a)=(a+a')/2, S(b)=(b+b')/2 with the quaternion variables to achieve the same thing. Here we go. We first multiply our equation on the left by a: a(ax+xb) = ac aax + axb = ac (1) We then multiply our equation on the right by b': (ax+xb)b' = cb' axb' + xbb' = cb' (2) Adding (1)+(2): aax + axb + axb' + xbb' = ac + cb' aax + ax(b+b') + bb'x = ac + cb' aax + (b+b')ax + bb'x = ac + cb' (aa + (b+b')a + bb') x = ac + cb' x = (aa + (b+b')a + bb')^-1 (ac + cb') which makes sense, so long as the coefficient of x isn't zero. But the equation is symmetrical, so we can work instead with a' & b: a'(ax + xb) = a'c a'ax + a'xb = a'c (3) (ax + xb)b = cb axb + xbb = cb (4) Adding (3)+(4): a'ax + a'xb + axb + xbb = a'c + cb a'ax + (a'+a)xb + xbb = a'c + cb xaa' + x(a+a')b + xbb = a'c + cb x (aa' + (a+a')b + bb) = a'c + cb x = (a'c + cb) / (bb + (a+a')b + aa') We thus have solutions for our Sylvester equation. But wait, there's more! I claim that the case where N(a)=N(b) is a most important case. Suppose that we have a 'particular' solution xp of ax+xb=c, i.e., a xp + xp b = c. If there are any other solutions of ax+xb=c, then they would involve 'homogeneous' solutions xh of the 'homogeneous' equation ax+xb=0, i.e., a xh + xh b = c. We can now add these two equations: a xp + xp b + a xh + xh b = c a(xp+xh) + (xp+xh)b = c which shows that xp+xh is also a solution of our original equation ax+xb=c. Let's better understand homogeneous solutions xh: a xh + xh b = 0 <=> a xh = - xh b => N(a xh) = N(- xh b) <=> N(a)N(xh) = N(-1)N(xh)N(b) <=> N(a)N(xh) = N(xh)(b) Cancelling N(xh) from both sides, we find that we require N(a)=N(b) for nontrivial solutions xh. But if N(a)=aa'=N(b)=bb'=N, then our particular solutions from above look like: x = (aa + (b+b')a + bb')^-1 (ac + cb') = (aa + (b+b')a + a'a)^-1 (ac + cb') = ((a + (b+b') + a')a)^-1 (ac + cb') = ((a+a'+b+b')a)^-1 (ac + cb') = ((2 S(a) + 2 S(b))a)^-1 (ac + cb') = (1/a) (ac + cb') / (2S(a)+2S(b)) = ((1/a)ac + (1/a)cb') / (2S(a)+2S(b)) = (c + a'cb'/N(a)) / (S(a)+S(b))/2 = (c + a'cb'/N) / (2S(a)+2S(b)) And x = (a'c + cb) / (bb + (a+a')b + aa') = (a'c + cb) / (bb + (a+a')b + b'b) = (a'c + cb) / ((b + (a+a') + b')b) = (a'c + cb) / ((a+a'+b+b')b) = (a'c + cb) / ((2S(a)+2S(b))b) = (a'c + cb) / (b(2S(a)+2S(b)) = (a'c/b + cb/b) / (2S(a)+2S(b)) = (a'cb'/bb' + c) / (2S(a)+2S(b)) = (a'cb'/N + c) / (2S(a)+2S(b)) = (c + a'cb'/N) / (2S(a)+2S(b)) In other words, *both particular solutions are identical* ! If we wish, we can write our equation in a slightly more symmetrical fashion for understanding eigenvectors: (S(a)+S(b)) x = (c + a'cb'/N)/2 Here, the coefficient of the unknown quaternion x is simply the scalar number (S(a)+S(b)); the 4D quaternion (c+a'cb'/N)/2 gives the appropriate *direction* in 4D. This form also makes obvious the problem that occurs when S(a)=-S(b); we can't divide by zero, so our 'solution' formula doesn't work in this case. We utilize (c+a'cb'/N)/2 instead of c+a'cb'/N because we will soon find that the map c -> (c+a'cb'/N)/2 is a useful *projection* operation, which preserves quaternions c which are oriented in certain directions.