Solving for Ax=B where x>=0

Open discussion regarding features, bugs, issues, vendors, etc.

Solving for Ax=B where x>=0

Postby Farfalo » Sat Oct 13, 2012 11:26 pm

Hi guys,

I have been looking at LAPACK for a few days now, and I have it working to an extent, however I can't seem to find a function that does exactly what I need. I need to solve Ax=B, where I have a 6x6 matrix A with one right-hand side vector and keep x>=0. I was able to use DGESV in order to solve for x, but the simplest solution in most cases involves some negative numbers, which in my case causes some issues as it involves chemical additions to water (you can't add negative amounts! :lol:).

I was wondering if there is a function in LAPACK or a way to set DGESV up differently so that x is greater than or equal to 0?

Thank you
Farfalo
 
Posts: 1
Joined: Sat Oct 13, 2012 11:15 pm

Re: Solving for Ax=B where x>=0

Postby fletchjp » Tue Oct 16, 2012 6:44 am

I think the problem is most likely somewhere in your data.

Try checking a hand solution of the problem.

John
fletchjp
 
Posts: 8
Joined: Wed Oct 03, 2012 10:13 am

Re: Solving for Ax=B where x>=0

Postby CyLith » Wed Oct 17, 2012 2:34 pm

If your A matrix is square and full rank, then there is only a single unique solution, and if the solution you get from Lapack has negative entries, then as John said, the problem is with your data being inconsistent with expectations. If your A matrix is rank deficient (has a null space), then there are an infinite number of solutions, and you will have to select one by some additional criteria (like positivity). You can find the least norm solution using the built in Lapack routines, and also find the nullspace of A with a QR decomposition; together these two things let you parameterize the entire space of solutions.
CyLith
 
Posts: 37
Joined: Sun Feb 08, 2009 7:23 am
Location: Stanford, CA

Re: Solving for Ax=B where x>=0

Postby Julien Langou » Wed Oct 17, 2012 2:53 pm

Hello,

Yep, there is no such routine in LAPACK. You want to do something called "nonnegative least squares".

The reference algorithm is given in Lawson, C.L. and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, 1974, Chapter 23, p. 161.

Here is a more recent code for it: http://hesperia.gsfc.nasa.gov/~schmahl/nnls/nnls.for

And I did not know about this but just found that there is a CUDA / Open MP implementation for the code:
http://www.cs.umd.edu/~yluo1/Projects/NNLS.html
See: Yuancheng Luo and Ramani Duraiswami, "Efficient Parallel Non-Negative Least Squares on Multicore Architectures", SIAM Journal on Scientific Computing 2011.

Cheers,
Julien.
Julien Langou
 
Posts: 734
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Solving for Ax=B where x>=0

Postby serkanulu » Mon Feb 11, 2013 10:26 pm

This is useful, "nonnegative least squares" (http://hesperia.gsfc.nasa.gov/~schmahl/nnls/nnls.for). i firstly understood, "least squares" and i thinking, imaginary values for example (ax + bi). i^2 (i square) is equal to -1 (minus one). But it's my fault, you defined (x >= 0). i think sleepy.
serkanulu
 
Posts: 1
Joined: Mon Feb 11, 2013 10:16 pm
Location: University of Izmir


Return to User Discussion

Who is online

Users browsing this forum: Bing [Bot], Yahoo [Bot] and 1 guest

cron