Create and solve different combinations of linear equation systems using data stored in several matrices in python

Question

I want to give an easy example for what I would like to do. I have some data which is calculated prior to this step and stored in matrices. In this easy case let’s take only two 2x2 matrices. They could look like this:

A1=np.array([[1, 2], [0.5, 1.5]])


A2=np.array([[0.5, 1.2], [1.3, 2]])

Also I have the solution vector b, that could look like:

b=[4, 3]

Now I would like to solve all possible combinations of linear equation systems Ax=b that can be created from the data stored in the matrices A1 and A2. In both matrices the data of the columns are connected to each other. So the possible linear equations systems to be solved from my example would look like:

LES1:

1 * x1 + 0.5 * x2 = 4

0.5 * x1 + 1.3 * x2 = 3

LES2:

1 * x1 + 1.2 * x2 = 4

0.5 * x1 + 2 * x2 = 3

LES3:

2 * x1 + 0.5 * x2 = 4

1.5 * x1 + 1.3 * x2 = 3

LES4:

2 * x1 + 1.2 * x2 = 4

1.5 * x1 + 2 * x2 = 3

With using np.linalg.solve(Ax, b) I get four sets of x1 and x2 values which will be evaluated in a next step to find the best combination of the input data matching a certain requirement. I know I can solve this problem using nested for loops to slice the input arrays A1 and A2 to create the linear equations systems LES1 to LES4 and solve them one after the other.

But in my real application I need to handle much more data than in this example which means I have more than two matrices and they are also bigger which will create much more possible combinations of linear equation systems that need to be solved.

I would be very very happy if someone could give me a hint to solve this problem without using for loops to hop through the columns of the matrices. Or maybe there is a possible way to use for loops in an efficient manner?

Thanks to all of you in advance

David

Update:

Mihails answer solves the problem I have described perfectly. However, when using my “real life data” two problems occurred which I have unfortunately not considered in my question and the described example.

1. The shape of the matrices is not always the same. Although the number of rows is identical in all matrices, the number of columns vary. As for example:

A1=np.array([[4, 2, 3, 5, 4], [1, 0.5, 5, 3, 1], [2, 3, 5, 4, 2]])
A2=np.array([[4, 2, 4], [8, 4, 1], [2, 8, 9]])
A3=np.array([[4, 8, 3, 2], [5, 6, 4, 5], [1, 2, 5, 7]])

Still I want to find all combinations of linear equation systems like mentioned in my actual question. Is there any way to adapt the given solution to this task? I tried to solve this problem reading the itertools documentation, but I couldn’t figure it out since it is not possible to create an array of arrays with this data like you have done in your answer when creating array A. And combining them in a List will give me a problem when it comes to the iteration with i.

2. After proceeding the data I need to know which input dataset fit best to my requirements. Means which combination of columns give the best result. Since the values in the columns are calculated in an earlier step from some measurement data I would like to carry on the name of the measurement data they are based on (Can be extracted from the Data in an earlier step). I know I can give column names when using pd.DataFrame for storage, but I am not sure if they can be used to perform the task I asked for. Thank you again and sorry for the obviously to basic “Mini example”

Answer 1

An elegant and simple approach is using itertools.product to avoid hardcoding loops. It takes multiple lists and computes the cartesian product. Down below you can see how it can compute the combinations of an arbitrary number of columns.

In terms of performance, zip and product are using generators, so they won't fill your memory directly but generate each element on the fly.

On the other hand, if you need to compute really huge amounts of data you could think about distributing the load to multiple cores/machines.

import numpy as np
from itertools import product

A1=np.array([[1, 2], [0.5, 1.5]])
A2=np.array([[0.5, 1.2], [1.3, 2]])
A = np.array([A1, A2])
b=[4, 3]

combine_ith_column = lambda A, i: product(*A[:, i])
combine_all_columns = lambda A: [combine_ith_column(A, i) for i in range(A.shape[1])]
A_space = zip(*combine_all_columns(A))

for Ax in A_space:
    print (np.linalg.solve(Ax, b))