




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>, >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.tudarmstadt.de (Peter
Spellucci) wrote: > 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
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
(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:
> 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
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.
(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:
> 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
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:
> On Jun 6, 2:18 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 > > (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 > 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
I think I may have found my problem. dgelsd requires a 1D matrix for the A matrix, not a 2D matrix.
(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:
> 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
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:
> 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
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...
(snip)
> 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...
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:
> byaa @yahoo.com wrote: > (snip) > > 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... > 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
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:
> On Jun 6, 5:26 pm, glen herrmannsfeldt <g @ugcs.caltech.edu> wrote: > > byaa@yahoo.com wrote: > > (snip) > > > 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... > > 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 > 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); > } Hide quoted text  >  Hide quoted text  >
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:
> On Jun 6, 6:30 pm, byaa @yahoo.com wrote: > > On Jun 6, 5:26 pm, glen herrmannsfeldt <g@ugcs.caltech.edu> wrote: > > > byaa@yahoo.com wrote: > > > (snip) > > > > 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... > > > 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 > > 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); > > } Hide quoted text  > >  Hide quoted text  > > > 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
Greg, Dave, thanks for your time. That did solve my problem.
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
On Jun 8, 5:21 am, byaa@yahoo.com wrote:
> 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
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:
> On Jun 8, 5:21 am, byaa @yahoo.com wrote: > > 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 > 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.
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>, >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





