Hi Ed,
On Fri, Jun 11, 2021 at 06:30:51PM +0000, Ed . wrote:
> 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?
I think the problem we are discussing is just a matter of
convention. Thus, I insist that it is not very important how this is
solved, as long as it is consistent.
As I understand it, Linear Algebra is related to the
linear transformation of vectors. The transformations are represented
by matrices and the vectors by (linear) arrays. Then a transformation T
acting on a vector V is another vector W=T V which is represented by the
product of the matrix T_{ij} by the array V_j, i.e., W_i=T_{ij} V_j
(summing over the repeated index j), where T_{ij} is the element on
the i-th row and j-th column. This would correspond to the usual inner
product of PDL, i.e., $W=$T->inner($V). However, vectors could also be
represented by one-column matrices. Then the product would be given by
the ordinary matrix product, not the inner product, i.e., $Wc=$T x
$Vc. The advantage of this approach is that you represent matrices and
vectors alike. Furthermore, if you build a matrix $Vmc from many
columns, you can transform all columns simultaneously through the
matrix product and obtain a matrix $Wmc out of the resulting columns,
i.e., $Wmc=$T x $Vmc.
Thus, currently I believe that the ideal solution would be for lu_backsub
to mimic the behavior of the matrix product 'x'. If you feed it a 1D
array it should complain and refuse to work. This shouldn't be
troublesome as one can convert from vector to column matrix by adding a
dummy dimension. Then the solution of a system of equations would also
be a column matrix. If the right hand side of the equation has n
columns, it should correspond to n systems of equations with the same
coefficients. Then the answer would also have n columns, one for each
solution. This way, the expected matrix, lu decomposition, right hand
side and solutions would all be 2D matrices. Then, one could do
ordinary threading over any extra dimensions in the matrix of
coefficients, their lu decompositions and/or right the hand sides of
the system of equations.
pdl> $m=sequence(2,2)
pdl> $x=sequence(1,2)
pdl> p $x
=>[[0],[1]] #edited
pdl> p $m x $x
=>[[1],[3]]
pdl> p $m x $x->dummy(2,2)
=>[[[1],[3]],[[1],[3]]]
pdl> p $m->dummy(2,2) x $x
=>[[[1],[3]],[[1],[3]]]
pdl> p $m->dummy(2,2) x $x->dummy(2,2)
=>[[[1],[3]],[[1],[3]]]
Actually, that is the current behavior of lu_backsub, except that it
represents a set of vectors as a matrix with a row per vector instead
of a column per vector for both the input and output. Thus, the need
for transpositions.
pdl> p lu_backsub(lu_decomp($m), $x->transpose)->transpose
=>[[0.5],[0]] #edited
pdl> p lu_backsub(lu_decomp($m), $x->transpose->dummy(2))->transpose
=>[[[0.5],[0]]]
pdl> p lu_backsub(lu_decomp($m), $x->transpose->dummy(2,2))->transpose
=>[[[0.5],[0]],[[0.5],[0]]]
pdl> p lu_backsub(lu_decomp($m->dummy(2,2)), $x->transpose)->transpose
=>[[[0.5],[0]],[[0.5],[0]]]
pdl> p lu_backsub(lu_decomp($m->dummy(2,2)),
$x->transpose->dummy(2,2))->transpose
=>[[[0.5],[0]],[[0.5],[0]]]
The problem with my original suggestion that 1D inputs should yield 1D
outputs would be that it would lead to an ambiguity in the
interpretation of a 2D or higher dimensional input. In
pdl> $m=random(3,3,4);
pdl> $y=random(3,4);
pdl> $x=lu_backsub(lu_decomp($m), $y)->transpose
should $x be interpreted as 4 1D vectors $x(:,($j)) $j=0..3
each of which obeys the system
$m(:,:,($j)) x $x(:,$j)->transpose=$y(:,$j)->transpose
(i.e., four systems of equations, each with a different 3x3 coefficients,
different 3-vector right hand sides and different 3-vector solutions for
each $j) or should it be interperted as a set of matrices $x(:,:,($j))
each of which obeys
$m(:,:,($j)) x $x(:,:,($j))->transpose=$y(:,:)->transpose
that is, four system of equations, each with different coefficients
matrix but with the same 3x4 matricial right hand
side and with four 3x4 matricial solutions.
By demanding that the input be 2D at least there would be no ambiguity
Thus my suggestion would be to keep the current behavior of lu_backsub
except that it should die if the input is a 1D vector, or complain, or
the documentation should say that 1D vectors would be converted to 2D
matrices. The documentation should also say that the second dimension
of the right hand side is not a threading dimension.
The signature should then look like
lu(m,m); perm(m); b(m,n); [o] out(m,n)
with an explicit n for the right hand side and the output.
Sorry for the long and somewhat obscure answer.
Regards,
Luis
> 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 :-)
(See 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
>
--
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 | moc...@fis.unam.mx /\_/\__/
GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB
_______________________________________________
pdl-general mailing list
pdl-general@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/pdl-general
