Python: решение матричного уравнения Ax = b, где A содержит переменные

Question

Я надеюсь решить матричное уравнение вида Ax = b, где

A = [[a,0,0],
     [0,a,0],
     [0,0,a],
     [1,1,0],
     [0,0,1]]

x = [b,c,d]

b = [0.4,0.4,0.2,0.1,0.05]

Если бы я знал значение a, я сделал бы это с помощью модуля numpy.linalg.lstsq. Точно так же, если бы я думал, что решение было точным, я бы использовал симпози решатель. Тем не менее, решение должно соответствовать наименьшему квадрату для некоторых данных, и решение для симпли обычно говорит мне, что решения не существует.

В основном я ищу версию numpy.linalg.lstsq, которая позволяет мне объявлять одну или несколько переменных в A. Есть ли готовый модуль, который занимается этими типами проблем?

Answer 1

Я придумал что-то, что вы можете использовать для решения общего случая. Также вам не нужно указывать границы или диапазоны.

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

"Решить нелинейную задачу наименьших квадратов с границами для переменных."

#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

Answer 2

Кажется, что проблема, которую вы опубликовали, имеет точное аналитическое решение:

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

Суммирование первых 2 уравнений и деление их на 4-е приводит к:

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

=> a=8

Таким образом, вы можете, если вы знаете вектор b, вы всегда можете найти для ab, ac ad, а затем использовать 4-й для решения для a.

Answer 3

Как вы упомянули в комментариях, проблема слишком ограничена, поэтому нет точного решения.

С другой стороны, проблема с использованием np.linalg.lstsq напрямую заключается в том, что у вас есть переменная в матрице A.

Один из возможных обходных путей - определить диапазон для переменной a и найти решение для параметра c в зависимости от него, что я имею в виду:

import numpy as np
import pandas as pd

# Define the list of results
R = []

b = np.array([0.4, 0.4, 0.2, 0.1, 0.5])

# Here I picked a (-10,10) range for a with a step of 0.01 
for a in np.arange(-10, 10, 0.01):
    # Define A given a
    A = np.array([[a,0,0], [0,a,0], [0,0,a], [1,1,0], [0,0,1]])
    # Solve with least-squares and store both a and the result
    R.append((a, *np.linalg.lstsq(A,b)[0]))

# Convert solutions to a dataframe
results = pd.DataFrame(R, columns=list('abcd'))

results
          a         b         c         d
0    -10.00 -0.038235 -0.038235 -0.014851
1     -9.99 -0.038271 -0.038271 -0.014861
2     -9.98 -0.038307 -0.038307 -0.014871
3     -9.97 -0.038343 -0.038343 -0.014880
4     -9.96 -0.038379 -0.038379 -0.014890
    ...       ...       ...       ...
1995   9.95  0.040395  0.040395  0.024899
1996   9.96  0.040355  0.040355  0.024870
1997   9.97  0.040315  0.040315  0.024840
1998   9.98  0.040275  0.040275  0.024811
1999   9.99  0.040236  0.040236  0.024782

Оказавшись в этой форме, вы можете увидеть, как переменные b, c и d изменяются вместе с a, нанося на график результаты:

results.set_index('a').plot(figsize=(16,8))

Обратите внимание, что b==c для любого a, и поэтому вы не видите ни одного следа на графике.

Answer 4

Вот уравнения, которые вы хотите решить:

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⎦

Мы можем получить решение наименьших квадратов, используя псевдообратную форму:

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                          ⎦

Answer 5

Поскольку проблема нелинейная (включая такие продукты, как ab a c и ad), вам понадобится какой-то нелинейный метод. Вот один итерационный метод. в каждой итерации мы будем оценивать b, c, d, используя линейный метод наименьших квадратов, а затем используем 3 оставшихся уравнения для оценки a. и мы сделаем это до сближения.

Обозначим переменную, которую мы пытаемся оценить, как x1, x2 и x3, поскольку b - это вектор с правой стороны. используя эту запись, мы имеем

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

, что означает:

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

, а остальные уравнения представляют собой простую линейную систему:

A_d * x = b[2:]

Матрица A_d имеет вид аналогично матрице A без первых 3 строк, т.е. A_d = A [2 :,:].

Теперь мы можем использовать итеративное решение:

# 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   

Вы можете выберите один из многих других критериев сходимости.