Solving overdetermined system in numpy when the value of one variable is already known
Question: I'm trying to solve an overdetermined system in Python, using the
`numpy.solve` function. I know the value of one of the variables and I know
that in theory I can find a unique solution for the system if I can somehow
plug in that known value.

My system is of the form `AxC=B`. The variables are split into two groups, one
group of N variables and one of T variables (although this does not matter for
the math). A is a `(T*N x T+N)` matrix, C is the variables vector, of length
`(T+N)`, and B is a vector of length `(T*N)`.

How do I tell `numpy.solve` (or another function in Python, but please don't
recommend least squares, I need the unique, exact solution, which I know
exists) to use the known value of one of the variables?

A simple example of my system would be:

    
    
    |1 0 0 1 0|     |n1|     |B1|
    |1 0 0 0 1|     |n2|     |B2|
    |0 1 0 1 0|  X  |n3|  =  |B3|
    |0 1 0 0 1|     |t1|     |B4|
    |0 0 1 1 0|     |t2|     |B5|
    |0 0 1 0 1|              |B6|  
    

The values of the elements of B would of course be known, as well as the value
of one of the variables, let's say I know that `t1=1`. The dots don't mean
anything I just put them there so the characters wouldn't bunch up.


Answer: As @Foon pointed out, the canonical way to do this is to subtract a column.

However, on a side note, as your problem is overdetermined, you _have_ to use
a method such as least squares. By definition, if it's an overdetermined
problem, there is no "unique, exact solution". (Otherwise it would be even-
determined - A square matrix.)

That aside, here's how you'd go about it:

Let's take your example equation:

    
    
    |1 0 0 1 0|     |n1|     |B1|
    |1 0 0 0 1|     |n2|     |B2|
    |0 1 0 1 0|  X  |n3|  =  |B3|
    |0 1 0 0 1|     |t1|     |B4|
    |0 0 1 1 0|     |t2|     |B5|
    |0 0 1 0 1|              |B6|  
    

As you noted, this is overdetermined. If we know one of our "model" variables
(let's say `n1` in this case), it will be even more overdetermined. It's not a
problem, but it means we'll need to use least squares, and there isn't a
completely unique solution.

So, let's say we know what `n1` should be.

In that case, we'd re-state the problem by subtracting `n1` multiplied by the
first column in the solution matrix from our vector of observations (This is
what @Foon suggested):

    
    
    |0 0 1 0|     |n2|     |B1 - n1|
    |0 0 0 1|     |n3|     |B2 - n1|
    |1 0 1 0|  X  |t1|  =  |B3 - 0 |
    |1 0 0 1|     |t2|     |B4 - 0 |
    |0 1 1 0|              |B5 - 0 |
    |0 1 0 1|              |B6 - 0 | 
    

Let's use a more concrete example in numpy terms. Let's solve the equation `y
= Ax^2 + Bx + C`. To start with, let's generate our data and "true" model
parameters:

    
    
    import numpy as np
    
    # Randomly generate two of our model variables
    a, c = np.random.rand(2)
    b = 1
    x = np.linspace(0, 2, 6)
    
    y = a * x**2 + b * x + c
    noise = np.random.normal(0, 0.1, y.size)
    y += noise
    

First, we'll solve it _without) the knowledge that `B = 1`. We could use
`np.polyfit` for this, but to lead into the next bit, we'll use a lower-level
approach:

    
    
    # I'm a geophysist, so I tend to use Gm=d instead of Ax=b
    G = np.column_stack([x**2, x, np.ones_like(x)])
    d = y
    
    m, _, _, _ = np.linalg.lstsq(G, d)
    
    print "Ideally, this would be 1: ", m[1]
    

As you can see, we'll get something close to, but not quite 1. In this case (I
didn't set the seed, so this will vary), the model parameters returned are

    
    
    [ 0.13392633,  0.97217035,  0.33645734]
    

While the true parameters are:

    
    
    [ 0.14592752,  1.        ,  0.31349185]
    

Now let's take the fact that we know `b` exactly into account. We'll make a
new `G` with one less column and subtract that column times `b` from our
observations (`d`/`y`):

    
    
    G = np.column_stack([x**2, np.ones_like(x)])
    d = y - b * x
    m, _, _, _ = np.linalg.lstsq(G, d)
    

Now `m` is `[a, c]` and we've solved for those two variables using our
knowledge of `b`.

