Thread Subject: underdetermined linear equation system

hello,

I have a big an very urgent problem. I have to solve an equation system in this way:

I have 50 unknows, and 8-15 equations. I have an initiative solution, so I'm looking for an iterative method to solve the system.

Since I'm not a mathematician, I googled a lot for an appropiate method, but I couldn't find nothing.
I tried the mldivide function of the matlab, but it is not good for me. the solutions must be all positive, and matlab gives me a lot of zeros (I know that's the way it solves the system with mldivide).

So I'm stuck.
can anybody help me please?

bets regards
arnold

"arnold szilagyi" <john.doe.nospam@mathworks.com> wrote in message <hasu5b$di4$1@fred.mathworks.com>...
> hello,
>
> I have a big an very urgent problem. I have to solve an equation system in this way:
>
> I have 50 unknows, and 8-15 equations. I have an initiative solution, so I'm looking for an iterative method to solve the system.
>
> Since I'm not a mathematician, I googled a lot for an appropiate method, but I couldn't find nothing.
> I tried the mldivide function of the matlab, but it is not good for me. the solutions must be all positive, and matlab gives me a lot of zeros (I know that's the way it solves the system with mldivide).
>

What you must specify how you want your solution to be, among the infinity number of solutions of

A*x = b % A is (8 x 50) matrix

For example, if you tell
|x| is minimal, % where |x| is sqrt(sum(x.^2))

then the solution is
x = pinv(A)*b

If you want solution to be all positive, then you need to add

A*x = b
x >= epsilon

where epsilon is a small positive number. The positive constraints alone are *not* enough solve the system. You might want

min J(x) := |x|^2 such that
A*x = b
x >= epsilon

This is quadratic programming. If you have optilmization toolbox, the function QUADPROG will do that for you.

Bruno

To complement my previous posts, the two other popular constraints are:

(P1)
min L1(x) := sum |x(i)|
A*x = b
x >= epsilon

(P2)
min Linf(x) := max{|x(i)|}
A*x = b
x >= epsilon

Both problems (P1) and (P2) can be solved with "linear programming" using LINPROG (also in optimization toolbox).

Bruno

"arnold szilagyi" <john.doe.nospam@mathworks.com> wrote in message <hasu5b$di4$1@fred.mathworks.com>...
> hello,
>
> I have a big an very urgent problem. I have to solve an equation system in this way:
>
> I have 50 unknows, and 8-15 equations. I have an initiative solution, so I'm looking for an iterative method to solve the system.
>
> Since I'm not a mathematician, I googled a lot for an appropiate method, but I couldn't find nothing.
> I tried the mldivide function of the matlab, but it is not good for me. the solutions must be all positive, and matlab gives me a lot of zeros (I know that's the way it solves the system with mldivide).
>
> So I'm stuck.
> can anybody help me please?
>
> bets regards
> arnold

There are infinitely many solutions to your problem.
So MATLAB chooses an arbitrary one. If you have
additional information (i.e., that the solution must
be non-negative) then you must supply that piece
of information!

Use a tool that will allow you to do this. One such
tool is linprog. A second is lsqnonneg. I'd expect to
see different solutions from the two, but since there
are infinitely many solutions, what else can you
expect?

Note that if you use linprog, you will need to provide
a linear objective function.

John

thank you for your quick and precious reply.

first of all I tried the quadprog function with the min J(x) := |x|^2 cost function.
but all that matlab gives to me is zeros.
here is the code that doesn't work for me (maybe I'm missing something):

A=[0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0;0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1;1 0 0 1 1 1 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0;0 0 1 0 0 0 1 1 1 1 0 1 0 0 0 1 0 0 0 0 1 1 0 1 0 1 1 1 1 1 0 0 0 0 0 0 1 0 1 0 1 1 0 0 1 1 1 0 0 0;1 1 1 1 0 1 1 1 0 0 1 1 1 1 1 0 0 0 1 1 0 0 0 1 0 1 1 0 1 1 1 1 1 1 1 1 0 0 0 1 0 0 1 0 0 1 1 1 0 0;0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 1 1 0 0 1 1 1 0 1 0 0 1 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0 1 1 0 0 0 1 1;0 0 1 0 1 1 1 0 0 0 0 1 0 1 1 1 1 0 1 1 1 0 1 1 0 1 0 1 1 1 0 0 1 1 0 1 1 0 1 0 0 1 0 0 0 1 0 1 1 0;1 1 0 1 0 0 0 1 1 1 1 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 0 0 0 1 1 0 0 1 0 0 1 0 1 1 0 1 1 1 0 1 0 0 1;1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 1 0
0;0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0;0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;0 1 1 1 0 0 0 1 1 0 0 1 1 1 0 0 0 0 0 1 1 0 0 1 0 1 0 0 1 1 1 1 1 0 1 0 1 1 0 0 1 0 1 0 0 1 1 0 1 0;0 0 0 0 0 0 1 0 1 1 1 0 0 1 0 1 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 1 0 1 0 1 1 1 0 0 0 0 0 1 0;1 1 0 1 1 0 0 1 0 0 0 1 1 0 1 0 0 1 1 0 0 1 0 0 1 1 1 1 0 0 1 1 1 1 1 0 0 1 0 1 0 0 0 0 1 1 1 1 0 1;0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0];

B=[312;351;396;506;930;580;891;621;585;360;84;575;496;756;288];

n=51;
H=toeplitz([1 zeros(1,n-2)],[1 zeros(1,n-2)]);

f=zeros(50,1);

lb=zeros(50,1);

ub=ones(3,1);
ub=ub*100;

[x,fval,exitflag,output,lambda] = quadprog(H,f,A,B,[],[],lb,ub)


this is the situation: I look for 50 numbers which's average satisfies some constraints.
for example:
(a1+a2+a4)/3=20
(a1+a3)/2=60

so I made the A matrix corresponding to the groups of numbers and I made B in this way (adapted to my little example):

A=[1 1 0 1; 1 0 1 0];
B=[20*3; 60*2];

I made a diagonal matrix for H with the values of 1 on the diagonal, and put 0 to the f, because I just use the x^2 in the cost function.
I added lower and upper limits and called the function.
the answer is 50 zero.
I don't know what am I missing.


after this I tried the linprog function
every matrix is the same and I call

[x,fval,exitflag,output,lambda] = linprog(f,A,B,[],[],lb,ub)

now it gives me nice numbers as a solution. but they don't satisfy the average conditions.
Maybe made a mistake to give to F zeros. but I don't know what should I try.

any suggestions?

regards
arnold

"arnold szilagyi" <john.doe.nospam@mathworks.com> wrote in message <hata91$ra1$1@fred.mathworks.com>...
maybe I'm missing something

> [x,fval,exitflag,output,lamda] = quadprog(H,f,A,B,[],[],lb,ub)
 
Yes READ careful the doc of quadprog/linprog, you mess inputs between inequalities and equality constraints.

Bruno

> Yes READ careful the doc of quadprog/linprog, you mess inputs between inequalities and equality constraints.
>
> Bruno

quote:
"x = quadprog(H,f,A,b,Aeq,beq) solves the problem above while additionally satisfying the equality constraints Aeq*x = beq."

this means that my A and B matrices are actually Aeq an beq???

if so what is the right syntax, because I can not figure it out.

arnold

>
> this means that my A and B matrices are actually Aeq an beq???
>
> if so what is the right syntax, because I can not figure it out.
>

Read even *more* carefully the doc.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hatmob$rm9$1@fred.mathworks.com>...
> >
> > this means that my A and B matrices are actually Aeq an beq???
> >
> > if so what is the right syntax, because I can not figure it out.
> >
>
> Read even *more* carefully the doc.
>
> Bruno

I'm sorry but I obviously miss something.
I read a lot of times the doc, and I made the following conclusion:

I have to use in this way:
[x,fval,exitflag,output,lamda] = quadprog(H,f,[],[],A,B,lb,ub)

but I tried and it gives me zeros.

Please help me.

regards
arnold

The model you have defined is infeasible, ie. not underdetermined but overdetermined and thus infeasible in quadprog

 output.message

ans =

The equality constraints are dependent.
The system of equality constraints is not consistent.
The inequality constraints may or may not be satisfied.
There is no feasible solution.


Just try to solve the linear systems and you'll see it directly

>> x = A\B;
>> A*x-B

ans =

   -7.3750
   -7.3750
   -7.3750
   -7.3750
   12.7500
   12.7500
   11.7500
   11.7500
  -17.1250
  -17.1250
  -17.1250
  -17.1250
         0
   -0.0000
   -0.0000

"arnold szilagyi" <john.doe.nospam@mathworks.com> wrote in message <hav1qg$mh7$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hatmob$rm9$1@fred.mathworks.com>...
> > >
> > > this means that my A and B matrices are actually Aeq an beq???
> > >
> > > if so what is the right syntax, because I can not figure it out.
> > >
> >
> > Read even *more* carefully the doc.
> >
> > Bruno
>
> I'm sorry but I obviously miss something.
> I read a lot of times the doc, and I made the following conclusion:
>
> I have to use in this way:
> [x,fval,exitflag,output,lamda] = quadprog(H,f,[],[],A,B,lb,ub)
>
> but I tried and it gives me zeros.
>
> Please help me.
>
> regards
> arnold

"arnold szilagyi" <john.doe.nospam@mathworks.com> wrote in message <hav1qg$mh7$1@fred.mathworks.com>...


>
> but I tried and it gives me zeros.

Your A matrix is rank 12 and B is not in Image A. So there is never a x such that A*x=B. Garbage in garbage out.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hav5cs$beq$1@fred.mathworks.com>...
> "arnold szilagyi" <john.doe.nospam@mathworks.com> wrote in message <hav1qg$mh7$1@fred.mathworks.com>...
>
>
> >
> > but I tried and it gives me zeros.
>
> Your A matrix is rank 12 and B is not in Image A. So there is never a x such that A*x=B. Garbage in garbage out.
>
> Bruno

can you explain me this this, please?

I can not work it out, even if I reduce the rank of A?

thank you
arnold

"arnold szilagyi" <john.doe.nospam@mathworks.com> wrote in message <havpas$8c2$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hav5cs$beq$1@fred.mathworks.com>...
> > "arnold szilagyi" <john.doe.nospam@mathworks.com> wrote in message <hav1qg$mh7$1@fred.mathworks.com>...
> >
> >
> > >
> > > but I tried and it gives me zeros.
> >
> > Your A matrix is rank 12 and B is not in Image A. So there is never a x such that A*x=B. Garbage in garbage out.
> >
> > Bruno
>
> can you explain me this this, please?
>
> I can not work it out, even if I reduce the rank of A?
>

Your case is similar to this simple example:
A=[0 0 0;
     1 0 0]
B=[1; 2];

In this example rank(A) is 1. And span(A) = { A*x: x in R^3 } = { (0,a): a in R }. I let you check this statement.

Since B(1)=1 in the above example, B is not in span(A) so there is no x such that A*x=B. Therefore when you feed the constraint A*x = B into QUADPROG and LINPROG, it cannot find any x. The "feasible set" is empty. Thus there is no solution (you should check the output flag, which is described in the doc) returned by QUADPROG.

In short, your system is not only underdetermined, but it's also rank degenerated.

Now back to *your* problem with 15 x 50 matrix size, do this to check the rank of A and alseo to see if you encountered similar problem:

>> rank(A)

ans =

    12

>> B

B =

   312
   351
   396
   506
   930
   580
   891
   621
   585
   360
    84
   575
   496
   756
   288

>> BPrj=A*pinv(A)*B % Projection of B on span(A)

BPrj =

  304.6250
  343.6250
  388.6250
  498.6250
  942.7500
  592.7500
  902.7500
  632.7500
  567.8750
  342.8750
   66.8750
  557.8750
  496.0000
  756.0000
  288.0000

>> Err = B-BPrj % Err is the projection of B on orthogonal(span(A))

Err =

    7.3750
    7.3750
    7.3750
    7.3750
  -12.7500
  -12.7500
  -11.7500
  -11.7500
   17.1250
   17.1250
   17.1250
   17.1250
    0.0000
    0.0000
    0.0000

>> norm(Err)/norm(B) % <- should be close to 0 for B in span(A)

ans =

    0.0205

% Bruno