Approximation to unsolvable system of equations
Let
$$x_1 := \log_2 (a) \qquad \qquad x_2 := \log_2 (b) \qquad \qquad x_3 := \log_2 (c) \qquad \qquad x_4 := \log_2 (d)$$
Thus, the original system of bilinear equations can be transformed into a system of linear equations
$$\begin{bmatrix} 1 & 0 & 1 & 0\\ 1 & 0 & 0 & 1\\ 0 & 1 & 1 & 0\\ 0 & 1 & 0 & 1\end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3\\ x_4\end{bmatrix} = \begin{bmatrix} -1\\ \log_2 \left(\frac 5 6\right)\\ -3\\ -1\end{bmatrix}$$
Find the least-squares solution $(\hat x_1, \hat x_2, \hat x_3, \hat x_4)$ and then compute the corresponding $(\hat a,\hat b,\hat c,\hat d)$.
In MATLAB:
>> A = [1 0 1 0; 1 0 0 1; 0 1 1 0; 0 1 0 1]
A =
1 0 1 0
1 0 0 1
0 1 1 0
0 1 0 1
>> b = [-1; log2(5/6); -3; -1]
b =
-1.0000
-0.2630
-3.0000
-1.0000
Unfortunately, the matrix is singular:
>> det(A)
ans =
0
We then find one solution to the "normal equations":
>> x = (A' * A)\(A' * b)
x =
1.7194
0.3509
-3.0351
-1.6667
The corresponding $(\hat a,\hat b,\hat c,\hat d)$ is, thus:
>> y = 2.^x
y =
3.2930
1.2754
0.1220
0.3150
We introduce a new matrix:
>> Y = y * y'
Y =
10.8437 4.1997 0.4017 1.0372
4.1997 1.6266 0.1556 0.4017
0.4017 0.1556 0.0149 0.0384
1.0372 0.4017 0.0384 0.0992
The absolute error is:
>> Y(1:2,3:4) - [1/2,5/6;1/8,1/2]
ans =
-0.0983 0.2039
0.0306 -0.0983
which is quite large.
Hint: If you stuff the variables into vectors, your problem becomes:
$$M = \left[\begin{array}{c}A\\B\end{array}\right]\left[\begin{array}{cc}C&D\end{array}\right] = \left[\begin{array}{rr} \frac{1}{2} & \frac 56\\\frac 1 8&\frac 1 2 \end{array}\right]$$
This means we want to find the best rank 1 match to the matrix.
Thanks for @Ian s comment which clarified nicely that you can use for example the Singular Value Decomposition ( SVD ) to find the best match.
The SVD says $M = U\Sigma V^*$ where the values in the diagonal $\Sigma$ are the singular values. If we pick the column in U and the row in $V^*$ which correspond to the largest singular value, we can identify the $A,B,C$ and $D$ above.
Some Matlab/Octave code to do the SVD approximation:
M = [1/2,5/6;1/8,1/2];
[U,S,V]=svd(M);
U(:,1)*S(1,1)*V(:,1)' % rank 1 SVD appr. as index 1 contains largest s.v.
$$\left[\begin{array}{cc} 0.445505865263938&0.861513345600555\\ 0.230380036476643&0.445505865263938 \end{array}\right]$$
abs(M-U(:,1)*S(1,1)*V(:,1)') %and the absolute error
$$\left[\begin{array}{cc} 0.0544941347360618&0.028180012267222\\ 0.105380036476643&0.0544941347360618 \end{array}\right]$$
We can see the error distributed so that each element has a rather different error. It differs $10.5/2.8 \approx 3.74$ times smallest compared to largest. As mentioned in comments, if we are not OK with this distribution of the error, we may want to try minimize the error according to some other norm.
Another way you could try to approach this problem by using a least squares estimator.
Assume $$F(A,B,C,D)=(AC-1/2)^2+(AD-5/6)^2+(BC-1/5)^2+(BD-1/2)^2$$
Now calculate the gradient of $F$ $$\nabla F=\left[\dfrac{\partial F}{\partial A},\cdots,\dfrac{\partial F}{\partial D}\right]$$
Then set the gradient to zero. Solve the resulting nonlinear system using Newton-Raphson algorithm (you should also check the positiveness of the Hessian).