# Solve a non-trivial matrix equation for an unknown matrix

I am solving the matrix equation to pass from the Fokker-Planck to the Langevin equation. Specifically: $$$$D=BB^T$$$$ where $$D= \left(\begin{matrix} l x_1(1-r_1) + e_1 r_1 & 0 & l x_1(1-r_1) & 0 \\ 0 & l x_2(1-r_2) + e_2 r_2 & 0 & l x_2(1-r_2) \\ l x_1(1-r_1) & 0 & l x_1+2D_{11}x_1 + D_{12}x_1 +D_{21}x_2 c_{11}r_1 +c_{12}r_2 & D_{12}x_1 +D_{21}x_2 \\ 0 & l x_2(1-r_2) & D_{12}x_1 +D_{21}x_2 & l x_2+2D_{22}x_2 + D_{12}x_1 +D_{21}x_2 c_{21}r_1 +c_{22}r_2 \end{matrix}\right)$$

where $$B$$ is symmetric. I am trying to solve the equation for $$\mathbb{B}$$ as follows:

Dmatrix = {{e1 r1 + l x1 (1 - r1), 0, l x1 (1 - r1), 0}, {0,      e2 r2 + l x2 (1 - r2), 0, l x2 (1 - r2)}, {l x1 (1 - r1), 0,      l x1 + 2 D11 x1 + c11 r1 + D12 x1 + D21 x2 + c12 r2,      D12 x1 + D21 x2}, {0, l x2 (1 - r2), D12 x1 + D21 x2,      l x2 + D21 x2 + 2 D22 x2 + D12 x1 + c21 r1 + c22 r2}};

B = Array[Y, {4, 4}];
For[i = 1, i < 5, i++, For[j = i, j < 5, j++, B[[j, i]] = B[[i, j]]]];

sol = Solve[Dmatrix == B . Transpose[B], Flatten[B]]; B /. Flatten[sol] // MatrixForm


But I am not finding a proper solution. Mathematically I know that more than a solution will be possible, since these are 16 quadratic equations. Can anyone help me?

• This is a polynomial system of 16 equations in 16 unknowns with several parameters. I think no chance to find its symbolic solution. Sep 19, 2023 at 16:12
• Finding (what is effectively) the square root of a matrix, is a notoriously difficult problem. One may use a particular decomposition, eg Cholesky etc (see the help). Sep 19, 2023 at 17:07

Doing a Cholesky decomposition, gives the answer practically instantaneously:

Dmatrix = {{e1 r1 + l x1 (1 - r1), 0, l x1 (1 - r1), 0}, {0,
e2 r2 + l x2 (1 - r2), 0, l x2 (1 - r2)}, {l x1 (1 - r1), 0, l x1 + 2 D11 x1 + c11 r1 + D12 x1 + D21 x2 + c12 r2,D12 x1 + D21 x2}, {0, l x2 (1 - r2), D12 x1 + D21 x2,l x2 + D21 x2 + 2 D22 x2 + D12 x1 + c21 r1 + c22 r2}};

B = CholeskyDecomposition[Dmatrix] //. Conjugate[X_] -> X // Simplify

(* {{Sqrt[e1 r1 + l x1 - l r1 x1], 0, -((l (-1 + r1)
x1)/Sqrt[e1 r1 - l (-1 + r1) x1]), 0}, {0, Sqrt[
e2 r2 + l x2 - l r2 x2], 0, -((l (-1 + r2) x2)/Sqrt[e2 r2 -
l (-1 + r2) x2])}, {0, 0, Sqrt[ c11 r1 + c12 r2 + 2 D11 x1
+ D12 x1 + l x1 + ( l^2 (-1 + r1)^2 x1^2)/(-e1 r1 + l (-1 +
r1) x1) + D21 x2], ( D12 x1 + D21 x2)/Sqrt[ c11 r1 + c12 r2
+ 2 D11 x1 + D12 x1 + l x1 + ( l^2 (-1 + r1)^2 x1^2)/(-e1
r1 + l (-1 + r1) x1) + D21 x2]}, {0, 0, 0, \[Sqrt](c21 r1 +
c22 r2 + D12 x1 + D21 x2 + 2 D22 x2 + l x2 - (D12 x1 + D21
x2)^2/( c11 r1 + c12 r2 + 2 D11 x1 + D12 x1 + l x1 + ( l^2
(-1 + r1)^2 x1^2)/(-e1 r1 + l (-1 + r1) x1) + D21 x2) + (
l^2 (-1 + r2)^2 x2^2)/(-e2 r2 + l (-1 + r2) x2))}} *)


The result looks a bit ugly, but it works:

Transpose[B].B - Dmatrix // Simplify

(* {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}} *)


In the above I have made the very bold assumption that everything is real (note the Conjugate[X_] -> X part, so caveat emptor!

Finally, as mentioned previously this decomposition is not unique, however is extremely fast (0.00081 secs as reported by RepeatedTiming on my machine).

This is only half the answer. I am able to solve for D, but to get B from D seems another problem.

Note, B is symmetric. But then B.Transpose is B.B. And this is also symmetric.

We may solve Dmatrix == D, where D is symmetric:

D = Array[Y, {4, 4}];
For[i = 1, i < 5, i++, For[j = i, j < 5, j++, D[[j, i]] = D[[i, j]]]];

sol = Solve[Dmatrix == B, DeleteDuplicates[Flatten[B]]][[1]];
B2 = B /. sol

{{e1 r1 + l x1 - l r1 x1, 0, -l (-1 + r1) x1, 0}, {0,
e2 r2 + l x2 - l r2 x2, 0, -l (-1 + r2) x2}, {-l (-1 + r1) x1, 0,
c11 r1 + c12 r2 + 2 D11 x1 + D12 x1 + l x1 + D21 x2,
D12 x1 + D21 x2}, {0, -l (-1 + r2) x2, D12 x1 + D21 x2,
c21 r1 + c22 r2 + D12 x1 + D21 x2 + 2 D22 x2 + l x2}}


But now we need to solve B2 == B.B. For this MMA takes longer than I am willing to wait.

• After a while B=MatrixPower[B2, 1/2]  evaluates... Sep 19, 2023 at 17:02
• @UlrichNeumann: Can you kindly present the result and timing? I ask it because that command crashes my comp in short time. Sep 19, 2023 at 18:43
• @user64494 Evaluation B=MatrixPower[B2, 1/2]; takes 25s on Mathematica v12.2, but I can't display the huge result in finite time Sep 20, 2023 at 8:45
• @user64494 B//LeafCount gives 182264749, probably to much to display without crash? Sep 20, 2023 at 9:15