Forcing conditions in a solution of an under-determined system

[Telemachus]

Hi there. I am solving an under-determined system of equations. The solver I've written seems to work fine, but I would like to have some control on the solution I obtain. From the infinite set of solutions of my system, I would like to use only those that fulfill some conditions. One of those conditions is that the elements of the solution should ##x_k\geq-1##. The other condition I would like it to meet is more difficult to explain, but basically it has to do with that there should be some relation between some sets of elements in the solution. Basically, the element ##x_1## should be equal to the element, let's say ##x_M##, and ##x_{M-c}## where c and M are constants and so on.

Is there a way to meet these conditions? I am programming this in fortran, and I am using lapack to solve the under-determined system (I'm using the dgels subroutine). The program finds a solution, but sometimes the conditions are not fulfilled, and the solutions I obtain are physically nonsensical.

Bye there.

 

[StoneTemplePython]

Ok, there are lots of ways to solve an under-determined system of equations. Two common and popular ones involve coming up with a solution vector ##\mathbf x## in your equation of ##A\mathbf x = \mathbf b##, such that ##\mathbf x## has a minimum "length" -- using a 2 norm or a 1 norm.

The 2 norm approach is done via typical linear algebra solvers -- but they don't accommodate the constraints you're asking about directly. (There might be an indirect workaround though it seems tedious at best and it would depend on how detailed your constraints are)

On the other hand, the minimizing 1 norm approach is done via Linear Programming and is quite flexible. Supposing that the constraints you want are feasible (as in they don't preclude a feasible / valid solution), a linear programming approach can easily accommodate what you're asking for.

Long story short: consider setting this up as a linear program.

- - - -
note: I use Julia for Linear Programming -- I've never done it in Fortran.

It looks a bit dreary, but maybe something like this would be helpful for Fortran specifics:
http://www.ccom.ucsd.edu/~peg/papers/sqdoc.pdf

edit:
for avoidance of doubt: encoding a constraint that certain parameters in your solution vector should be equal can be done via typical matrix methods and solvers. If you have a nice LP modeling setup like in Julia's JuMP, I think it's nicer to do it there, but there is nothing preventing you from tacking on new rows to the bottom of your ##A## to enforce these constraints.

The real issue comes in dealing with inequality constraints like ##x_k \geq -1##. Linear Programming in some sense was designed to deal with things like this. 'Regular' matrix solvers are not. I would not recommend, in effect, creating Simplex from scratch via coming up with solutions via 'regular' solvers and doing lots of pivotting on your own... it is much better to use an off the shelf, highly efficient LP solver.

 

[Telemachus]

Ok, thanks a lot.