On Mon, Jun 5, 2017 at 1:21 PM, Craig DeForest <[email protected]>
wrote:
> I agree that there’s an issue here! It can easily become a sort of
> religious war like emacs vs vi. Many (most?) languages have chosen the
> route you suggest — the (row,column) mathematical style indexing is kept,
> and matrices are rendered backward. In the late Jurassic (was it Tuomas
> who wrote the ‘x’ operator?) PDL chose (column,row) to match (x,y). That
> helps folks who like to visualize matrices but hurts folks who prefer to
> think in indices. Folks who like to visualize matrices can see their 3x4
> PDL as a 3-column by 4-row matrix (what the mathematicians would call a
> 4x3). Folks who like to work with index notation have to remember that
> arrays are indexed in CR (XY) order, not RC (YX) order, which can be a pain
> to remember.
>
My proposal is that we compute with PDLs with matrix operations as if they
were matrices: (i,j) would correspond to (dim0,dim1) which is the same
ordering
used in PDL slice operations and the display operation would be *decoupled*
from the dimension ordering for computation (I've always found it annoying
in MATLAB when an image is displayed as a matrix rather than as an image).
Done this way matrix operations would use (dim0, dim1) ordering and image
operations would also use (dim0, dim1) ordering.
> Doing it the other way would have the converse problem: we’d be addressing
> matrices in the order chosen by the mathematics community, but
> visualizations be transposed. This includes basics like images, which are
> themselves a subclass of objects that might be considered “matrices” — and
> then the image processing community (like that guy DeForest) would
> complain, since images are universally addressed in XY order and not RC
> order.
>
I'm proposing that the data display *not* determine the computational use
of the PDL. This way we could keep natural operation ordering for linear
algebra while supporting equally the (x,y) view for image operations. Any
isomorphism between matrix and image representations can easily be
converted to an isomorphism between matrix and the image-transposed
representation as well.
>
> I suppose one could produce a PDL::Matrix subclass that transposes its
> first two dimensions, sidestepping the whole issue by producing a
> specialist object. One would have to think carefully about how that would
> interact with ordinary active dims and threading.
>
>
The current PDL::Matrix object forces computation on PDLs with row
major ordering to behave like column major ordering. If we chose this
split display/compute approach, PDL::Matrix would only need to have
a stringify/display that does column major rather than row major.
I think that would be much simpler and more robust than the current
dimension re-ordering behind the scenes.
Cheers,
Chris
> Best,
> Craig
>
>
> On Jun 5, 2017, at 11:07 AM, Chris Marshall <[email protected]>
> wrote:
>
> On Mon, Jun 5, 2017 at 11:29 AM, Craig DeForest <[email protected]
> > wrote:
>>
>>
>> I don't find PDL's matrix handling to be screwy. There's a
>> necessary wart between column-major addressing and row-major
>> addressing, sure but that's endemic to all languages. The
>> question is whether you want matrices to _naturally_render_
>> the way they would appear in mathematical notation, or whether
>> you want them to _naturally index_ the way they would appear
>> in mathematical notation. Some languages (e.g., IDL) chose the
>> latter; others (e.g., PDL) chose the former. You can't have both
>> without breaking the way arrays are rendered on-screen (so that,
>> e.g., "$a = pdl(1,2,3)" would render as a column vector, taking 5
>> lines of text.
>>
>>
> Here's an example of what I mean for screwy with matmult:
>
> pdl> $a34 = (10*random(3,4))->floor;
> pdl> p $a34
> pdl> p $a34->transpose
> pdl> $b42 = (10*random(4,2))->floor;
> pdl> p $b42
> pdl> p $a34 x $b42
> pdl> p $a34->transpose x $b42->transpose
> pdl> p +($a34->transpose x $b42->transpose)->transpose
> pdl> p $b42 x $a34
> pdl> p $b42
>
> which yields on pasting into a pdl2 shell session:
>
> pdl> $a34 = (10*random(3,4))->floor;
>
> pdl> p $a34
>
> [
> [8 8 5]
> [3 8 4]
> [4 6 2]
> [9 4 9]
> ]
>
>
> pdl> p $a34->transpose
>
> [
> [8 3 4 9]
> [8 8 6 4]
> [5 4 2 9]
> ]
>
>
> pdl> $b42 = (10*random(4,2))->floor;
>
> pdl> p $b42
>
> [
> [4 0 8 0]
> [9 7 9 6]
> ]
>
>
> pdl> p $a34 x $b42
> Runtime error: Dim mismatch in matmult of [3x4] x [4x2]: 3 != 2 at
> /cygdrive/c/Perl/local64/lib/perl5/cygwin-thread-multi/PDL/Primitive.pm
> line 265. PDL::matmult(PDL=SCALAR(0x6038d5dd8), PDL=SCALAR(0x6038d71b0),
> PDL=SCALAR(0x6038cb860)) called at /cygdrive/c/Perl/local64/lib/
> perl5/cygwin-thread-multi/PDL/Core.pm line 819
> PDL::__ANON__(PDL=SCALAR(0x6038d5dd8),
> PDL=SCALAR(0x6038d71b0), "") called at (eval 449) line 5
>
> pdl> p $a34->transpose x $b42->transpose
>
> [
> [ 64 183]
> [ 80 206]
> [ 36 145]
> ]
>
>
> pdl> p +($a34->transpose x $b42->transpose)->transpose
>
> [
> [ 64 80 36]
> [183 206 145]
> ]
>
>
> pdl> p $b42 x $a34
>
> [
> [ 64 80 36]
> [183 206 145]
> ]
>
>
> pdl> p $b42
>
> [
> [4 0 8 0]
> [9 7 9 6]
> ]
>
> My observation is that if you look at the PDL shapes and in
> memory dataordering, then $a34 is laid out *exactly* as a
> fortran a(3,4) array would be.
>
> Similarly, the $b42 is laid out *exactly* as the fortran b(4,2)
> array would be.
>
> The default display ordering for PDLs is essentially row
> major but an alternative way to view it might be as in
> memory order. PDL historically has used the display
> order to determine the axis order for matrix operations.
> The result is to multiply a [3,4] shape piddle by a [4,2]
> piddle one must either transpose all the arguments and
> then transpose the result, or reverse the order of the
> multiplicands.
>
> If instead we use the standard PDL dimension numbering
> reading from left to right we could have the matrix multiply
> operation defined as in mathematics or as tensor sums:
>
> a_ij b_jk => c_ik
>
> The only "problem" is that the default display is tranposed.
> If we were to define a matrix print as one that displays in
> column major, then we get good math indexes and a
> consistent display.
>
>
>
>> As it stands now, "print ($vec=pdl([1,5]))" yields "[1 5]",
>> which is nice for more general contexts than just matrix
>> operations. Also, "print ($m=pdl([1,1],[0,1])" yields
>>
>> [
>> [1 1]
>> [0 1]
>> ]
>>
>> which is also nice: items are rendered in the most natural way.
>> That choice forces row-major ordering, which is the opposite of
>> the convention the mathematics community chose. (Can't blame
>> `em even Ben Franklin screwed up the sign convention for
>> electric charge#) The wart is that, if you want to hit a
>> column vector with $m, you have to say "$m_vec = $m x ($vec->(*1))"
>> to explicitly make $vec a column vector the 0 dim works along
>> a row, not a column.
>>
>> I don't see a clean way around the dichotomy between natural
>> array rendering and use (which is row-major) and matrix notation
>> (which is column-major). A transpose has to happen somewhere
>> either in the way arrays are rendered (column-major specification,
>> which makes matrices work great but screws up other applications)
>> or in the way that they are created (row-major specification,
>> which makes column vectors more cumbersome to use but makes other
>> applications more convenient).
>>
>
> Yes, there will be a dichotomy. My proposal is that lets put the
> inconsistency in the display and not in the computation. If you
> ignore how we currently display PDL data, there is *no* difference
> between a fortran array dimensions and indexing (m,n) and PDL
> dimensions and indexing [m,n] so why force the opposite convention?
>
> Am I the only one who finds it tricky to correctly convert
> matrix-vector algebra equations into the correct PDL ops?
> If the difference were moved to the print from the weird ordering
> of the matmult routine, then all that would be needed for matrix
> "sanity" would be to change the stringify operation for PDLs
> used as matrices or vectors.
>
> Cheers,
> Chris
>
>
>>
>> Sorry if this rambles, I wrote this before dashing off to a meeting.
>>
>> Best,
>> Craig
>>
>>
>> On Jun 4, 2017, at 3:36 PM, Chris Marshall <[email protected]>
>> wrote:
>>
>> But...
>>
>> If you look at the dimension ordering in slicing
>> where dim(0) is the left-most index, then PDL is
>> actually using the same memory ordering as with
>> fortan: dim0 iterates first, then dim1 increases,
>> then dim2....
>>
>> In fact sequence(3,4) is in memory in this
>> order: (0,0), (1,0), (2,0), (0,1), (1,1), (2,1),
>> (0,2), (1,2), (2,2), (0,3), (1,3), (2,3) which
>> is exactly the memory order of a(3,4) in fortran.
>>
>> It seems the issue is that *displaying* the
>> data uses C ordering. If we were to display
>> the data as if transposed, then PDL would seem
>> to me to be a column major storage system.
>>
>> I wonder what would happen if matrix operations
>> used the natural dimension order rather than
>> the imposed C ordering? It would get rid of
>> all the nasty transposes in the matrix multiplication
>> and things the tensor sums would compose naturally.
>>
>> Am I the only one who thinks PDL for matrix ops is
>> a bit screwy---but for no good reason?
>>
>> --Chris
>>
>> On 6/4/2017 16:27, Grégory Vanuxem wrote:
>>
>> Hi here,
>>
>> https://docs.oracle.com/cd/E19957-01/805-4940/z400091044d0/index.html
>>
>> Just for information.
>>
>> Now, an other thing to know.
>> Imagine I have a 2x2 matrix.
>>
>> If I write in the computer memory 4 integer’s (1-2-3-4) in one time, in C
>> this will be in matrix representation :
>>
>> 1 2
>> 3 4
>>
>> But in Fortran, like in mathematics :
>>
>> 1 3
>> 2 4
>>
>> So operations on these two representations are completely differents.
>>
>> Generaly computations on matrices are done on a very low level and use
>> directly the memory areas (no aware of indexing scheme).
>>
>> Hope that helps
>> __
>> Greg
>>
>> *De : *Luis Mochan <[email protected]>
>> *Envoyé le :*jeudi 1 juin 2017 20:06
>> *À : *[email protected]
>> *Objet :*Re: [Pdl-general] SVD
>>
>> Still confusing:
>> On Wed, May 31, 2017 at 03:36:33PM +1000, Karl Glazebrook wrote:
>> > column-major all the way down, as image processing came before matrix
>> ops and i in A(i,j) is naturally the x-axis.
>> but the x axis displays horizontally, as row, and when matrices are
>> multiplied i is interpreted as column index, i.e., as index along row.
>> > i.e. A[0,1] is followed by A[1,0] in memory
>> You mean A[0,0] is followed by A[1,0] in memory, right?
>> So memory is arranged as in fortran arrays (first index fastest), but
>> the interpretation as row and column indices is different.
>> Regards,
>> Luis
>>
>>
>>
>> >
>> >
>> > Karl
>> >
>> >
>> >
>> > > On 22 May 2017, at 6:32 am, Chris Marshall <[email protected]>
>> <[email protected]> wrote:
>> > >
>> > > Please ignore the following. Just mark me confused...
>> > >
>> > > --Chris
>> > >
>> > > On 5/21/2017 15:51, Chris Marshall wrote:
>> > >> The row-major and col-major for PDL has always
>> > >> confused me since, AFAICT the PDL dimensions and
>> > >> slicing syntax actually are column major in
>> > >> memory but we print them out in row-major.
>> > >>
>> > >> Maybe one of the original PDL developers could
>> > >> give an explanation (of the history at least!).
>>
>>
>
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
pdl-general mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/pdl-general
