Python: Solving Matrix Equation Ax = b where A contains variables

These are the equations you want to solve:

In [39]: a, b, c, d = symbols('a, b, c, d', real=True)                                                                            

In [40]: A = Matrix([[a,0,0], 
    ...:      [0,a,0], 
    ...:      [0,0,a], 
    ...:      [1,1,0]]) 
    ...: x = Matrix([[b],[c],[d]]) 
    ...: b = Matrix([[0.4],[0.4],[0.2],[0.1]])                                                                                    

In [41]: A*x - b                                                                                                                  
Out[41]: 
⎡ a⋅b - 0.4 ⎤
⎢           ⎥
⎢ a⋅c - 0.4 ⎥
⎢           ⎥
⎢ a⋅d - 0.2 ⎥
⎢           ⎥
⎣b + c - 0.1⎦

We can get the least squares solution using the pseudoinverse:

In [53]: A.pinv() * b                                                                                                             
Out[53]: 
⎡      ⎛ 2    ⎞                   ⎛ 2    ⎞            ⎤
⎢0.4⋅a⋅⎝a  + 1⎠     0.4⋅a     0.1⋅⎝a  + 1⎠      0.1   ⎥
⎢────────────── - ───────── + ──────────── - ─────────⎥
⎢   4      2       4      2     4      2      4      2⎥
⎢  a  + 2⋅a       a  + 2⋅a     a  + 2⋅a      a  + 2⋅a ⎥
⎢                                                     ⎥
⎢      ⎛ 2    ⎞                   ⎛ 2    ⎞            ⎥
⎢0.4⋅a⋅⎝a  + 1⎠     0.4⋅a     0.1⋅⎝a  + 1⎠      0.1   ⎥
⎢────────────── - ───────── + ──────────── - ─────────⎥
⎢   4      2       4      2     4      2      4      2⎥
⎢  a  + 2⋅a       a  + 2⋅a     a  + 2⋅a      a  + 2⋅a ⎥
⎢                                                     ⎥
⎢                         0.2                         ⎥
⎢                         ───                         ⎥
⎣                          a                          ⎦

It seems like the problem you've posted has an exact analytical solution:

ab = 0.4,
ac = 0.4,
ad = 0.2,
b + c = 0.1

Summing the first 2 equations and dividing it by the 4th results in:

(1) + (2) => a(b+c) = 0.8
(4) b+c = 0.1

=> a=8

So you can if you know the vector b you can always solve for ab, ac ad and than use the 4th to solve for a.

I came up with something you can use for a general case solution. Also you dont have to specify bounds or ranges.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html

"Solve a nonlinear least-squares problem with bounds on the variables."

#You can rewrite matrix A as

[[1, 0, 0],
 [0, 1, 0],
 [0, 0, 1],   * a + 
 [0, 0, 0],
 [0, 0, 0]]

[[0, 0, 0],
 [0, 0, 0],
 [0, 0, 0],
 [1, 1, 0],
 [0, 0, 1]]


#Define your loss function
def lossf(M):
    Y = np.array([0.4,0.4,0.2,0.1,0.05])
    a,b,c,d = M
    A = np.array([[1,0,0],[0,1,0],[0,0,1],[0,0,0],[0,0,0]])*a
    A = A + np.array([[0,0,0],[0,0,0],[0,0,0],[1,1,0],[0,0,1]])
    y = A.dot(np.array([b,c,d]).T)
    return np.sum((Y - y)**2)

#An extra function to calculate the matrix in the end
def mat(M):
    Y = np.array([0.4,0.4,0.2,0.1,0.05])
    a,b,c,d = M
    A = np.array([[1,0,0],[0,1,0],[0,0,1],[0,0,0],[0,0,0]])*a
    A = A + np.array([[0,0,0],[0,0,0],[0,0,0],[1,1,0],[0,0,1]])
    y = A.dot(np.array([b,c,d]).T)
    return y

import numpy as np
from scipy.optimize import least_squares

#this takes some time ** does 100000 rounds which can be changed
#[0,0,0,0] = random initialization of parameters
result = least_squares(lossf,[0,0,0,0], verbose=1, max_nfev = 100000)


result.x #result for parameter estimate
[6.97838023, 0.05702932, 0.05702932, 0.02908934] # run code and 

mat(result.x)
#The non-linear fit
[0.39797228, 0.39797228, 0.20299648, 0.11405864, 0.02908934]

#Orignal 
[0.4,0.4,0.2,0.1,0.05]


#Also results for other matrix
#This requires changing the loss function
#For a permanent solution you can write a flexible loss function
Y = np.array([0.4,0.4,0.2,0.1])
result = least_squares(lossf,[0,0,0,0], verbose=1, max_nfev = 10000)


result.x
[7.14833526, 0.0557289 , 0.0557289 , 0.02797854]

mat(result.x)
[0.39836889, 0.39836889, 0.2       , 0.11145781]
#The results are very close to analytical solution

Since the problem is non linear (include products like ab ac and ad) you'll need some sort of non-linear method. Here is one iterative method. in each iteration we'll evaluate b,c,d using linear least square, and then use the 3 remaining equation to estimate a. and we'll do it until convergence.

I'll denote variable we're trying to estimate as x1, x2 and x3 since b is the vector on the right hand side. using this notation we have

a * x1 = b[0]
a * x2 = b[1]
a * x3 = b[2]

which implies:

a = sum(b[0:3])/sum(x)  # x = [x1, x2, x3]

and the remaining equations are are a simple linear system:

A_d * x = b[2:]

The matrix A_d is the same as the matrix A without the first 3 rows, i.e. A_d = A[2:, :].

Now we can use an iterative solution:

# choose any initial values

prev = np.ones(4)  ## this includes `a` as the last element

rel_tolerance = 1e-4  # choose a number
while True:
    x_sol = np.linalg.lstsq(A_d, b[2:])
    a_sol = sum(b[0:3])/sum(x_sol)
    sol = np.append(x_sol, a_sol)  

    diff_x = np.abs(sol - prev)
    # equivalent to max relative difference but avoids 0 division
    if np.max(diff_x) < rel_tolerance * prev:  
        break   

You can choose one of many other convergence criteria.