Hi Luis, Thank you for the kind words! I think your points are of course valid. However, I’m hoping you can help me understand your expectations a bit better, so the tests, and then the code (and maybe the docs!), can be made to match them better.
In terms of dimensions, my understanding was that where there are vectors in linear algebra, they are actually (in row * column notation) a 1 * n matrix, because I thought linear algebra was all about matrix multiplication, being inherently x*n matrices to be multiplied by n*y matrices. If I’m wrong on that, then could you clarify? And more generally, could you give some trivial examples of your expectations of what would happen with 1-D, 2-D and 3-D inputs? I believe I understand your expectations to be 1-D -> 1-D, but I’d like to understand your response to the “really a matrix” idea above :-) To try to answer your other questions, in PDL a 3*3*1 matrix will behave identically to a 3*3 one due to threading. This is partly because both cases really are a 9-item area of memory, and partly because threading will treat a dimension of length 1 and a dummy dimension the same, and expand them to match if you tried to add say a 3*3*1 to a 3*3*2, or similarly a 3*3 to a 3*3*2. By the way, the above is the source of a problem that I had with the “native complex” adaptation for Photonic – I got seg-faults calling a PDL::LinearAlgebra function with a 1x1 LU decomp, but (due to my mistakes with slicing), an 11x1 “B”. Because PDL (I believe) treats a dimension of 1 as a dummy, it cheerfully thought I must really have passed in an 11x11 LU, made the dimensions of both the LU and B read as 11, then the LAPACK routine obviously overran its bounds. I will soon make this instead throw an exception when it’s on a “[phys]” parameter, because there should never be seg-faults. Best regards, Ed From: Luis Mochan<mailto:[email protected]> Sent: 08 June 2021 04:54 To: perldl<mailto:[email protected]> Subject: Re: [Pdl-general] [Pdl-devel] lu_backsub Hi Ed, On Tue, Jun 01, 2021 at 12:57:51AM +0000, Ed . wrote: > That is what I was aiming to achieve (and have captured in hopefully-thorough > tests). My thinking is that the linear algebra stuff is founded on matrices, > with at least 2 dimensions. Therefore, getting something back that has at > least 2 dimensions (and always at least 2) seems to me the least surprising > outcome (especially given your previous note was apparently prompted by > surprise). > > What do you (and others) think? I believe the most important thing is consistency. However, linear algebra is about both matrices and vectors. Vectors may be represented by matrices, which makes the current behavior reasonable, but they may be represented as 1D arrays. I would have expected a 1D input to return a 1D output. My surprise was that lu_backsub returned 1D or 2D's depending on the entries of the matrix (on the value of the permutations array), so it was not a consistent behavior. Now it seems to be consistent when solving a single system of equations, which is very good. Nevertheless, I haven't tested but I worry about the behavior for other cases. For example, if instead of feeding lu_decomp a single matrix, it is fed with an array of matrices (as an ndarray of dimension >= 3), what output should I expect if I feed lu_backsub a single 1D vector, a matrix with a single row, a matrix with several rows? For example, a 3x3 matrix and a 3 vector currently return a 3x1 one row matrix. A 3x3x1 matrix and a 3 vector would also return a 3x1 matrix or a 3x1x1 matrix? A 3x3x2 matrix and a 3 vector, or a 3x1 vector, or a 3x2 matrix with two rows? The documentation says that lu_backsub takes a matrix and a vector and that it threads as usual, which seems not to be the case currently. So, in summary, I think any decision would be appropriate if it is clear, consistent and consistent with the documentation. But, everything else being equal, I would prefer that a 2D matrix and 1D vector yield a 1D output+the usual threading rules for additional dimensions. Best regards, Luis ps. And thanks for the enormous effort you are doing. > > Best regards, > Ed > > From: Luis Mochan<mailto:[email protected]> > Sent: 01 June 2021 01:46 > To: perldl<mailto:[email protected]> > Subject: Re: [Pdl-general] [Pdl-devel] lu_backsub > > I tested the new version and it seems to be consistent. > > So the current behavior is that if you feed lu_backsub > with a 1D n-vector you consistently get out a 2D nX1 row > matrix-vector, the same as if the input vector got an additional dummy > dimension. Is this the expected behavior? > > Best regards, > Luis > > > On Mon, May 31, 2021 at 09:55:30PM +0000, Ed . wrote: > > Hi Luis, > > > > The difference you’re seeing is due to a sort of mirage; when you transpose > > both outputs (suitable for matrix-multiplication) you’ll get the correct > > results, which are a 1x2 column vector. > > > > However, the lu_backsub code was quite complicated and in ascertaining the > > above, and also making it work “inplace” (did you know that the current > > “copy” method when “inplace” just turns off “inplace” and returns the given > > ndarray?) with tests, I had to make some fixes and changes. Now the dims of > > the results will all be the same. > > > > As soon as I’m confident these changes haven’t broken any other stuff, I’ll > > be releasing them. Thanks for the report! > > > > Best regards, > > Ed > > > > From: Luis Mochan<mailto:[email protected]> > > Sent: 29 May 2021 20:28 > > To: perldl<mailto:[email protected]>; > > perldl<mailto:[email protected]> > > Subject: [Pdl-devel] lu_backsub > > > > I found a strange behavior in lu_backsub from PDL::MatrixOps: > > > > pdl> my $M=pdl([1,2],[3,4]); > > pdl> my $M1=pdl([3,4],[1,2]); # interchange two rows > > pdl> my $y=pdl(1,1); > > pdl> my $x=lu_backsub($M->copy->lu_decomp, $y); > > pdl> my $x1=lu_backsub($M1->copy->lu_decomp, $y); > > pdl> p $x->info; > > PDL: Double D [2,1] > > pdl> p $x1->info; > > PDL: Double D [2] > > > > In both cases I'm solving the same 2x2 system of equations using a 1D > > vector as the right hand side, but in the second case I interchanged > > rows, so I expected the same solution. In the first case I get back a 'row > > vector', i.e., a matrix with one row. However, in the second case I > > got a real 1D vector, not a matrix. > > > > My guess is that the problem is related to the code around line 1298 > > of matrixops.pd, which contains some warnings: > > > > if($nontrivial) { > > if($y->is_inplace) { > > $y .= $y->dummy(1,$y->dim(0))->index($perm->dummy(1,1))->sever; > > # TODO: check threading > > $out = $y; > > } else { > > $out = > > $y->dummy(1,$y->dim(0))->index($perm->dummy(1,1))->sever; # TODO: check > > threading > > } > > } else { > > # should check for more matrix dims to thread over > > # but ignore the issue for now > > $out = ($y->is_inplace ? $y : $y->copy); > > } > > print STDERR "lu_backsub: starting with \$out" . pdl($out->dims) . > > "\n" if $PDL::debug; > > > > > > Best regards, > > Luis > > > > > > -- > > > > o > > W. Luis Mochán, | tel:(52)(777)329-1734 /<(*) > > Instituto de Ciencias Físicas, UNAM | fax:(52)(777)317-5388 `>/ /\ > > Av. Universidad s/n CP 62210 | (*)/\/ \ > > Cuernavaca, Morelos, México | [email protected] /\_/\__/ > > GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB > > > > > > _______________________________________________ > > pdl-devel mailing list > > [email protected] > > https://lists.sourceforge.net/lists/listinfo/pdl-devel > > > > -- > > o > W. Luis Mochán, | tel:(52)(777)329-1734 /<(*) > Instituto de Ciencias Físicas, UNAM | fax:(52)(777)317-5388 `>/ /\ > Av. Universidad s/n CP 62210 | (*)/\/ \ > Cuernavaca, Morelos, México | [email protected] /\_/\__/ > GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB > > > _______________________________________________ > pdl-general mailing list > [email protected] > https://lists.sourceforge.net/lists/listinfo/pdl-general > -- o W. Luis Mochán, | tel:(52)(777)329-1734 /<(*) Instituto de Ciencias Físicas, UNAM | fax:(52)(777)317-5388 `>/ /\ Av. Universidad s/n CP 62210 | (*)/\/ \ Cuernavaca, Morelos, México | [email protected] /\_/\__/ GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB _______________________________________________ pdl-general mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/pdl-general
_______________________________________________ pdl-general mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/pdl-general
