import numpy as np

# input_values
a = np.array([0.123456, 0.234567, 0.345678])
#outut_values
b = np.array([0.876543, 0.765432, 0.654321])


# A 3x3 Matrix multiplication $X a = b$ is a set of 3 equations
#
#   x_11 * a_1 + x_12 * a_2 + x_13 * a_3 = b_1
#   x_21 * a_1 + x_22 * a_2 + x_23 * a_3 = b_2
#   x_31 * a_1 + x_32 * a_2 + x_33 * a_3 = b_3
#
# There we have 3*3 = 9 unknown, but only 3 equations.  This doesn't work.  We
# have to cheat.  We just set 6 of the 9 unkown to 1.  For example every x_ij
# but the fist column
#
#   x_11 * a_1 + a_2 + a_3 = b_1
#   x_21 * a_1 + a_2 + a_3 = b_2
#   x_31 * a_1 + a_2 + a_3 = b_3
#
#   x_11 = b_1 - (a_2 + a_3) / a_1
#   x_21 = b_2 - (a_2 + a_3) / a_1
#   x_31 = b_3 - (a_2 + a_3) / a_1
#

X = np.ones((3,3))
X[0, 0] = (b[0] - (a[1] + a[2])) / a[0]
X[1, 0] = (b[1] - (a[1] + a[2])) / a[0]
X[2, 0] = (b[2] - (a[1] + a[2])) / a[0]

assert np.all( (np.dot(X, a) - b) < 1e-15 )  # approximately the same


# DOWNSIDE: Doesn't work if a_1 is zero or almost zero.  But this is easy to
# fix be picking the index i where |a_i| is the the greatest.

i_max = np.argmax(np.abs(a))
i_other = list(set(range(3)).difference([i_max]))

X = np.ones((3,3))
X[0, i_max] = (b[0] - np.sum(a[i_other])) / a[i_max]
X[1, i_max] = (b[1] - np.sum(a[i_other])) / a[i_max]
X[2, i_max] = (b[2] - np.sum(a[i_other])) / a[i_max]

assert np.all( (np.dot(X, a) - b) < 1e-15 )  # approximately the same