Home     |     .Net Programming    |     cSharp Home    |     Sql Server Home    |     Javascript / Client Side Development     |     Ajax Programming

Ruby on Rails Development     |     Perl Programming     |     C Programming Language     |     C++ Programming     |     IT Jobs

Python Programming Language     |     Laptop Suggestions?    |     TCL Scripting     |     Fortran Programming     |     Scheme Programming Language


 
 
Cervo Technologies
The Right Source to Outsource

MS Dynamics CRM 3.0

Fortran Programming Language

solving a system of linear equations


Hi,
In my experiments, I have a system of linear equations in 12
variables.  The 12 variables are essentially some context to my new
compression algorithm.  I have many such equations.

They are of the form

a1c1 + a2c2 + a3c3 ... + a12c12 = B1

I will typically have 500 of such equations.

My equations will not have an exact solution, but I am trying to find
the closest (best) solution for the variables a1, a2... a12, such that
the B vector is a close approximation of the input.

I am trying to figure out what function in LAPACK I can use to
accomplish this.  I beleive dgesv is incorrect since that will only
try to find the exact solution, which usually will not exist in my
data sets

B

In article <1181146539.625101.321@q19g2000prn.googlegroups.com>,

 byaa@yahoo.com writes:

 >Hi,
 >In my experiments, I have a system of linear equations in 12
 >variables.  The 12 variables are essentially some context to my new
 >compression algorithm.  I have many such equations.
 >
 >They are of the form
 >
 >a1c1 + a2c2 + a3c3 ... + a12c12 = B1
 >
 >I will typically have 500 of such equations.
 >
 >My equations will not have an exact solution, but I am trying to find
 >the closest (best) solution for the variables a1, a2... a12, such that
 >the B vector is a close approximation of the input.
 >
 >I am trying to figure out what function in LAPACK I can use to
 >accomplish this.  I beleive dgesv is incorrect since that will only
 >try to find the exact solution, which usually will not exist in my
 >data sets
 >
 >B
 >

dgelss  does exactly what you want
hth
peter

On Jun 6, 9:51 am, spellu@fb04373.mathematik.tu-darmstadt.de (Peter

I played with dgelss a little bit after your post.  I wasnt getting
the results I hoped for, so I inspected some of the output from
dgelss.

Even on simple equasions such as

a1c1 + a2C = C, where I have M such equasions (note that the constant
for the a2 variable is always the same as the RHS of the equasion),
dgelss would not find a proper solution, which in this case would be
simple a1 = 0.0 and a2 = 1.0

I must be missing something as to how this should work...

B

byaa@yahoo.com wrote:

(snip)

> Even on simple equasions such as
> a1c1 + a2C = C, where I have M such equasions (note that the constant
> for the a2 variable is always the same as the RHS of the equasion),
> dgelss would not find a proper solution, which in this case would be
> simple a1 = 0.0 and a2 = 1.0

It is usual to require the same number of equations as
unknowns for systems of linear equations.

In the case above, you have one equation and two unknowns.

Your simple solution requires that c1 not be zero, and also
the c1 not be C.

It might be that singular value decomposition would solve these
equations in the way you expect.  SVD allows for a different
number of equations from unknowns.

-- glen

On Jun 6, 1:55 pm, glen herrmannsfeldt <g@ugcs.caltech.edu> wrote:

I should have been clear.  In my case I have M equations with 2
unknowns.  What I was trying to say is that in each equation, the
coefficient of the a2 variable is always equal to the RHS in each of
the M equations.  I would have expected dgelsd to solve that properly.

B

In article <1181160062.294743.132@d30g2000prg.googlegroups.com>,

For M > 1,  a1 = 0, a2 = 1 will be a simultaneous solution of all the
equations, and very likely the only one.

byaa@yahoo.com wrote:

(snip)

>>>Even on simple equasions such as
>>>a1c1 + a2C = C, where I have M such equasions (note that the constant
>>>for the a2 variable is always the same as the RHS of the equasion),
>>>dgelss would not find a proper solution, which in this case would be
>>>simple a1 = 0.0 and a2 = 1.0

(snip)

> I should have been clear.  In my case I have M equations with 2
> unknowns.  What I was trying to say is that in each equation, the
> coefficient of the a2 variable is always equal to the RHS in each of
> the M equations.  I would have expected dgelsd to solve that properly.

And you have M=2?

If M is not 2 look at SVD.

If M is 2, is the coefficient matrix singular?  If it is,
even if the simple answer above satisfies the equations, you
won't get it from linear equation solvers.  Consider

3x+y=1
6x+2y=2

x=0, y=1 is a solution, but the matrix is singular and
linear system solvers won't give that solution.
Note that x=1, y=-2 is also a solution.

-- glen

On Jun 6, 2:18 pm, glen herrmannsfeldt <g@ugcs.caltech.edu> wrote:

Hi Glen, yes I understand that there are other solutions... and that
may be the case in my system.  However, the values for the variables
that dgelsd is returning does not correctly satisfy any of the input
equations.

Therefore my question is, since I know there exists atleast one best
solution that fully satisfies the equations, why is dgelsd still
chosing suboptimal answers?

B

On Jun 6, 1:42 pm, byaa@yahoo.com wrote:

I think I may have found my problem.  dgelsd requires a 1D matrix for
the A matrix, not a 2D matrix.

byaa@yahoo.com wrote:

(snip on using dgelsd)

It is SVD!  Ignore previous suggestions to use SVD.

> I think I may have found my problem.  dgelsd requires a 1D matrix for
> the A matrix, not a 2D matrix.

No, it can be 2D but the LDA argument must equal the first
dimension of A, especially if it isn't M.  You can do

       DOUBLE PRECISION A(100,100),B(101,11),S(100),WORK(10000)
       DOUBLE PRECISION RCOND
       INTEGER IWORK(10000)
       M=2
       N=2
       NRHS=1
       LDA=100
       LDB=101
       NWORK=10000
       A(1,1)=1.23
       A(1,2)=2.34
       A(2,1)=3.45
       A(2,2)=4.56
       B(1,1)=5.67
       B(2,1)=6.78
       RCOND=-1.
       CALL DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
      $                   WORK, LWORK, IWORK, INFO )

-- glen

On Jun 6, 3:23 pm, glen herrmannsfeldt <g@ugcs.caltech.edu> wrote:

Hmm... Why do you have NRHS = 1 but have a 2D matrix for B in your
example?
On Jun 6, 3:23 pm, glen herrmannsfeldt <g@ugcs.caltech.edu> wrote:

I am still having an issue... in my simplest case, I am doing the
following

In my equation AX = B,
my A matrix looks like this {0, -84, 0, -76}.  My B matrix looks like
this {-84 -76}

I call dgelsd as follows

one = 1;
rank = 0;
rcond = -1;
lwork = 10000;
info = 0;
M = 2;
N = 2;
dgelsd_(&M, &N, &one, A, &m, B, &N, S, &rcond, &rank, work, &lwork,
iwork, &info)

The resulting B matrix looks like this: {0.497506..., 0.450124688}

That seems wrong to me...

byaa@yahoo.com wrote:

(snip)

In Fortran arrays are stored with the first subscript varying fastest,
in C it is the other way around.  Your {0, -84, 0, -76} then corresponds
to the Fortran array elements A(1,1) A(2,1) A(1,2) A(2,2),
where M is the dimension of the leftmost subscript, and N the rightmost.

Does that help any?

If you use 2D array in C, be sure to use real 2D arrays, not
arrays of pointers to arrays.

double a[100][100],b[100][100];

where LDA will be the rightmost dimension in C.

-- glen

On Jun 6, 5:26 pm, glen herrmannsfeldt <g@ugcs.caltech.edu> wrote:

Maybe if you can point me how to rearrange my code, that will give me
an example... the example below is not working as expected:
#include <stdio.h>
#include <f2c.h>
#include <blaswrap.h>
#include <clapack.h>

static double A[100][100], B[100][100], S[10000], work[10000], rcond;
static long int M, N, one, rank, lwork, info, iwork[30000], lda, ldb;

int main()
{

        one = 1;
        rank = 0;
        rcond = -1;
        lwork = 10000;
        info = 0;
        lda = 100;
        ldb = 100;

        M = 2;
        N = 2;

        A[0][0] = 0;
        A[0][1] = -84;

        A[1][0] = 0;
        A[1][1] = -76;

        B[0][0] = -84;
        B[1][0] = -76;

    dgelsd_(&M, &N, &one, A, &lda, B, &ldb, S, &rcond, &rank, work,
&lwork, iwork, &info);

On Jun 6, 6:30 pm, byaa@yahoo.com wrote:

Changing these lines should do it...

        A[0][0] = 0;
        A[1][0] = -84;

        A[0][1] = 0;
        A[1][1] = -76;

        B[0][0] = -84;
        B[0][1] = -76;

Dave

On Jun 6, 5:48 pm, Dave Dodson <dave_and_da@Juno.com> wrote:

Greg, Dave, thanks for your time.  That did solve my problem.
On Jun 6, 1:14 pm, Virgil <vir@comcast.net> wrote:

My M = 2700 to be exact (that is, I have 2700 simultaneous equations
and exactly 2 variables, a1 and a2).  Furthermore, for each equation,
the coefficitent of a2 is always equal to the RHS value of that
equation.

So if I am not seeing a1 = 0 and a2 = 1, am I to conclude that I am
using dgelsd incorrectly?

B

On Jun 8, 5:21 am, byaa@yahoo.com wrote:

To avoid confusion let's call the variables by the names x1 and x2,
and the corresponding coefficients to be estimated by the names b1 and
b2.  For each element (observation) of M there is a corresponding
element of each of x1 and x2, and you seek to estimate b1 and b2, so
as (I presume) to minimize the sum of squared residuals ( the sum of
squared M - x1*b1 - x2*b2).  I searched the above conversation without
discovering any indication of what values are in the vectors x1 and
x2, and M; would you please clarify? It will be sufficient to say if
any of these vectors is constant across all observations, or varies,
and if it does vary, does it vary independently of the other vectors?
Also, are there variables included in this problem besides x1 and x2?

It could be that the results you have in hand are indeed the
solution.  This would be the case if M = some constant, x2 = vector of
1's, while x1 (and any other variables) varied across observations.
In that case, the least squares solution would set b2 = M's constant
value.

On Jun 8, 9:48 am, Mark Morss <mfmo@aep.com> wrote:

I should add that the least squares solution would also set b1 = 0.

b1 and b2 are >scalars<, of course.  Some of your posts seem to
suggest that you have different coefficients in different equations?!
I don't know what sort of analysis that would correspond to, but
certainly not to ordinary least squares estimation.

In article <1181162024.092009.150@i13g2000prf.googlegroups.com>,

 byaa@yahoo.com writes:

 >On Jun 6, 1:14 pm, Virgil <vir@comcast.net> wrote:
 >> In article <1181160062.294743.132@d30g2000prg.googlegroups.com>,
 >>
 >>
 >>
 >>  byaa@yahoo.com wrote:
 >> > On Jun 6, 1:55 pm, glen herrmannsfeldt <g@ugcs.caltech.edu> wrote:
 >> > > byaa@yahoo.com wrote:
 >>
 >> > > (snip)
 >>
 >> > > > Even on simple equasions such as
 >> > > > a1c1 + a2C = C, where I have M such equasions (note that the constant
 >> > > > for the a2 variable is always the same as the RHS of the equasion),
 >> > > > dgelss would not find a proper solution, which in this case would be
 >> > > > simple a1 = 0.0 and a2 = 1.0
 >>
 >> > > It is usual to require the same number of equations as
 >> > > unknowns for systems of linear equations.
 >>
 >> > > In the case above, you have one equation and two unknowns.
 >>
 > > > > Your simple solution requires that c1 not be zero, and also
 >> > > the c1 not be C.
 >>
 >> > > It might be that singular value decomposition would solve these
 >> > > equations in the way you expect.  SVD allows for a different
 >> > > number of equations from unknowns.
 >>
 >> > > -- glen
 >>
 >> > I should have been clear.  In my case I have M equations with 2
 >> > unknowns.  What I was trying to say is that in each equation, the
 >> > coefficient of the a2 variable is always equal to the RHS in each of
 >> > the M equations.  I would have expected dgelsd to solve that properly.
 >>
 >> > B
 >>
 >> For M > 1,  a1 = 0, a2 = 1 will be a simultaneous solution of all the
 >> equations, and very likely the only one.
 >
 >My M = 2700 to be exact (that is, I have 2700 simultaneous equations
 >and exactly 2 variables, a1 and a2).  Furthermore, for each equation,
 >the coefficitent of a2 is always equal to the RHS value of that
 >equation.
 >
 >So if I am not seeing a1 = 0 and a2 = 1, am I to conclude that I am
 >using dgelsd incorrectly?
 >
 >B
 >

 yes, necessarily. the svd computes under all possible solutions (in the
 case of rank deficinecy) the solution of shortest length.
 I saw somewhere in the code you declared the matrix as 100 times 100
 but had the leading dimension not =100.
 you must discern between actual dimension of your system (m and n)
 and the dimension you gave the matrix in their declaration _this is the LDA_

 your test system has the form

 [ c1(i), C ]   and the right hand side is always C

 if you have the vector of the c1(i) not constant, then this system will have
 full rank and you get a1=0 a2=C as the unique solution.
 if, as in your testcases, the matrix has first column all zeros,
 then of course a1 is arbitrary and a2=C, hence the solution  of minimal
 length gives a1=0, a2=C.
 and dgelss does exactly this.
 hth
 peter

Add to del.icio.us | Digg this | Stumble it | Powered by Megasolutions Inc