Is there a more legible way to solve sets of linear equations with numpy?

Question

I've got a set of 6 equations that I'd like numpy to solve for me. So I construct a 6x6 matrix of coefficients, and fill it in with various values. However, the code I end up writing to do this is quite illegible, and conveys little about the equations that I want to solve to the reader of my code.

For example, filling out the coefficients matrix looks something like this:

# Coefficients matrix
# Order of variables: w, X, Y, Z, s, t
A = np.mat( np.zeros((6,6)) )

A[0:3,0] = cam_inv[...,2]
A[0:3,1:4] = -np.identity(3)
A[3:6,1:4] = np.identity(3)
A[3:,4] = -eigvecs[...,0]
A[3:,5] = -eigvecs[...,1]

# Constants matrix (RHS of equation)
b = np.mat( np.zeros((6,1)) )
b[0:3,0] = -cam_inv[...,0:2] * point
b[3:,] = mean.T

res = np.linalg.solve(A,b)

(Where cam_inv, eigvecs, mean, and point are some other matrices computed elsewhere.)

Obviously the above code could have some more comments, but I feel that even with some comments, it'd still fail to really convey the underlying equations that are being solved. Is there a better way of feeding equations into the solver that results in code that is more legible?

Well it's entirely specific to your problem, isn't it? Is this some computer vision application? It might not lend itself to neat code but you could write out the elements of the matrix in a nice big comment above, with a reference to the equation/technique/paper you are implementing — YXD
Well indeed one could -- however, it feels like a problem that many people must have come across before, and so it seems like there should be a better interface to this situation. (And yes, it is a computer vision application.) — Rob Gilton

BKay BKay · Accepted Answer · 2013-04-05T13:06:14

The problem is that lines of A that represent the equalities don't have a one-to-one mapping to lines of code. What I do in my own work (Economics) is to have a function with a clear English name (or at least a single line of code, with no functional representation) for each of the rows of A. When necessary, I have a clear but slow or perhaps much longer version of the code that does the same thing as the code I ultimately use.

So for example (from Bretscher's Linear Algebra with Applications 1997, ex. 37 p. 29 a simple if unrealistic example), consider an economy with three industries, I1, I2, I3, each taking the other two industry outputs as inputs. What outputs should they produce to meet consumer and industrial demand?

A =np.zeros((3,3)) 
#Each unit of production by I1 requires 0.1 units of good 2 and .2 of good 3
A[:,0] = [0, 0.1, 0.2] 
#Each unit of production by I2 requires 0.2 units of good 1 and .5 of good 3
A[:,1] = [0.2, 0, 0.5] 
#Each unit of production by I3 requires 0.3 units of good 1 and .4 of good 2    
A[:,2] = [0.3, 0.4, 0] 
#The required production for consumers. 
b = np.array([320,150,90]).reshape(-1,1)
#The optimal production levels of x1, x2, and x3
res = np.linalg.solve(A,b)

It perhaps would be slower or less terse to do what I suggest, but it would be much clearer to read.

Is there a more legible way to solve sets of linear equations with numpy?

1 Answers