Beefy Boxes and Bandwidth Generously Provided by pair Networks
Just another Perl shrine

PDL::Complex question

by gmacfadden (Sexton)
on Mar 21, 2007 at 18:12 UTC ( #605900=perlquestion: print w/replies, xml ) Need Help??

gmacfadden has asked for the wisdom of the Perl Monks concerning the following question:

Manually multiplying
(1+i)(3+i) + (2+i)(2+i) = 5+8i
(1-2i)(3+i) +(2-i)(2+i) = 10-5i

These two equations can be expressed in terms of the matrix equation MX = B, where
/ 1+i 2+i \ / 3+i \ / 5+8i \ M = | | X = | | and B = | | \1-2i 2-i / \ 2+i / \ 10-5i /
One should therefore be able to use complex piddles to obtain B by multiplying M times X, or alternatively obtain X by multiplying the inverse of M times B using the following code:
#! /usr/bin/perl -w use warnings; use strict; use PDL; use PDL::Complex; my $matrixM = pdl [ [ 1+1*i, 2+1*i], [ 1-2*i, 2-1*i] ]; my $matrixB_assigned_value = pdl [ 5+8*i, 10-5*i ]; my $matrixX_assigned_value = pdl [ 3+1*i, 2+1*i ]; my ($matrixB_computed_value, $matrixX_computed_value); print "Via Perl Data Language\n"; print "List assigned_values of matrices:\n"; print "\$matrixM = ", $matrixM,"\n"; print "\$matrixX_assigned_value = ", $matrixX_assigned_value,"\n"; print "\$matrixB_assigned_value = ", $matrixB_assigned_value,"<br>\n"; print "List computed values of matrices:\n"; print "\$matrixB_computed_value = ", $matrixM x $matrixX_assigned_value,"\n"; print "\$matrixX_computed_value = ", $matrixM->inv x $matrixB_assigned_value,"\n"; exit(0);

When the above code is executed, the following output is created:
Via Perl Data Language List assigned_values of matrices: $matrixM = [ [ [1 1] [2 1] ] [ [ 1 -2] [ 2 -1] ] ] $matrixX_assigned_value = [ [3 1] [2 1] ] $matrixB_assigned_value = [ [ 5 8] [10 -5] ] <br> List computed values of matrices: $matrixB_computed_value = [ [ [5 2] [8 3] ] [ [-1 -1] [ 4 1] ] ] $matrixX_computed_value = [ [ [ 5 -13] [ 0 21] ] [ [ 5 -6] [ 0 -7] ] ]
Notice how the assigned values for M, X, and B dutifully print as expected. However the computed values for X and B are strange. Please can you explain what I need to change to correctly compute the values of X and B in my little example.

Replies are listed 'Best First'.
Re: PDL::Complex question
by drzowie (Initiate) on Apr 01, 2007 at 17:10 UTC
    The problem is that the 'x' operator is threading over the internal (real-valued) matrices rather than computing the correct complex products.

    That, in turn, is because you are using the pdl() constructor around your matrices. Thus they are not blessed into the PDL::Complex class, but only into the PDL class. The pdl() constructor is pretty permissive and turns just about anything into a PDL, but can be confused on occasion. In this case, pdl() recognizes a collection of nested list refs, some of which contain PDLs (or at least PDL subclasses). It is trying to do the Right Thing by putting your inner PDL dims in the deepest (fastest-running) indices of the final constructed PDL. That looks almost right, as PDL::Complex is implemented with an extra 0th dim that runs over {real|imaginary}.

    Because your matrices are PDLs and not PDL::Complexes, they are being multiplied differently than you expect: MatrixM is a 2x2x2, which is treated as a collection of two 2x2 matrices in (column, row, thread) order rather than a single complex matrix in ({r|i},column,row) order. MatrixX is a 2x2, which is being treated as a single matrix in (column, row) order rather than a complex row vector in ({r|i},column) order. You have a couple of ways around the complex issue:

    • (1) Use the correct constructor, "new PDL::Complex", rather than the incorrect one "pdl";
    • (2) Manually bless your PDLs back into PDL::Complex with 'bless($matrixM,"PDL::Complex");' (this is what PDL::Complex::new does anyway); or
    • (3) Assemble your matrices from independently constructed PDLs, as in '$matrix=pdl([5,4])+pdl([1,2])*i;'

    I think that you are about to hit another problem, which is that $matrixX is a row vector rather than (as you probably think, given that you put it on the right hand side of a matrix) a column vector.

    The row vector/column vector problem arises from the need to choose between (row,column) indexing order, which renders incorrectly on screen, or (column, row) indexing order, which is backwards but renders correctly on screen. PDL uses (column,row) order by default.

    To make a row vector, you say '$row = pdl(1,2,3)' or (if you prefer) '$row = pdl [1,2,3]'.

    To make a column vector you must insert a trivial 0th dim (of size 1), which you can do thusly:

    • (1) by nesting: $col = pdl([1],[2],[3])
    • (2) by transposition: $col = pdl(1,2,3)->transpose
    • (3) by dummy index: $col = pdl(1,2,3)->(*1)

    Incidentally, there's one more index pitfall to watch for (though it isn't yet in your example): sometimes you may want to thread across row vectors (for example, if you want to multiply a collection of N row vectors by a collection of N matrices that are stored in a CxRxN PDL). Then you should remember to insert the dummy dimension explicitly. For example, if your vectors are (1,2), (4,5), and (7,8) you should say:

    $matrices = sequence(2,2,3); $threadable_rows = pdl([1,2],[4,5],[7,8])->(:,*1,:); print $threadable_rows x $matrices;
    which will yield a different answer than
    $matrices = sequence(2,2,3); $wrong_rows = pdl([1,2],[4,5],[7,8]); print $wrong_rows x $matrices;
      Note from THE FUTURE: in PDL 2.040 in April 2021, "native complex" data (with types cfloat or cdouble) worked properly, eliminating the 0th "real or imaginary" dimension. That then made the x operator work in standard PDL with native complex numbers. Around the same time, PDL::LinearAlgebra and PDL::FFTW3 were made to also work with "native complex" data, allowing it to be used with those fast C/Fortran libraries.
Re: PDL::Complex question
by Zaxo (Archbishop) on Mar 22, 2007 at 06:53 UTC

    Column vectors should be piddled as,

    my $matrixB_assigned_value = pdl [ [5+8*i], [10-5*i] ]; my $matrixX_assigned_value = pdl [ [3+1*i], [2+1*i ] ];
    Untested, but I think that will give you the results you expect.

    After Compline,

      One should think so, but it doesn't.

      Conventional notation would have the 2x2 matrix left, the column vector to the right, so

      print " result: ", $matrixM x $matrixX_assigned_value;
      but that complains
      PDL: Dim mismatch in matmult of [2x2] x [2x1]: 2 != 1
      Swapping the factors
      print "should be: ", $matrixB_assigned_value; print " result: ", $matrixX_assigned_value x $matrixM;
      should be: [ [ [5 8] ] [ [10 -5] ] ] result: [ [ [5 4] ] [ [ 4 -5] ] ]
      The result looks like a column vector, but is numerically incorrect.

      BTW, it helps to cast the piddles into the PDL::Complex data type, using the cplx function. That overloads stringification so that the complex elements are printed in a x + i*y format instead of looking like pairs of reals. The numeric discrepancy remains, however.


Re: PDL::Complex question
by Anonymous Monk on Apr 05, 2007 at 06:29 UTC
    You can install PDL-LinearAlgebra if you want to do linear algebra over complexes. Here is an example (note that if you load PDL::LinearALgebra the x operator (matrix multiplication in PDL) is overriden by the BLAS matrix multiplication and thus is complex aware. Here is a perldl session:
    perldl> use PDL::Complex perldl> use PDL::LinearAlgebra perldl> $matrixM = cplx pdl [ [ 1+1*i, 2+1*i], [ 1-2*i, 2-1*i] ]; perldl> p $matrixM [ [1 +1i 2 +1i] [1 -2i 2 -1i] ] perldl> $matrixB_assigned_value = cplx pdl [ 5+8*i, 10-5*i ]; perldl> p $matrixB [ [ 5 +8i] [10 -5i] ] perldl> msolve($matrixM,$matrixB) [ [3+1i] [2+1i] ] [ [ 1 -2i 2 -1i] [-0.2 +0.6i 1.8 -0.4i] ] [2 2] 0 perldl> scalar msolve($matrixM,$matrixB) [ [3+1i] [2+1i] ] perldl> $matrixM->minv x $matrixB [ [3 +1i] [2 +1i] ]
    In the last result we see from the amount of space between the real and imaginary part that the result is less precise.
    ----------------------------------------------------------- perldl> help minv minv Compute inverse of a general square matrix using LU factorization. Supports inplace and threading. Uses getrf and getri or cgetrf and cgetri from Lapack and return "inverse, info" in array context. PDL(inv) = minv(PDL) my $a = random(10,10); my $inv = minv($a); Docs from /usr/lib/perl5/PDL/ perldl> help msolve msolve Solve linear system of equations using LU decomposition. A * X = B Returns X in scalar context else X, LU, pivot vector and info. B is overwrited by X if its inplace flag is set. Supports threading. Uses gesv or cgesv from Lapack. (PDL(X), (PDL(LU), PDL(pivot), PDL(info))) = msolve(PDL(A), PDL(B) +) my $a = random(5,5); my $b = random(10,5); my $X = msolve($a, $b); Docs from /usr/lib/perl5/PDL/

Log In?

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://605900]
Approved by andye
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others pondering the Monastery: (1)
As of 2022-10-01 02:37 GMT
Find Nodes?
    Voting Booth?
    I prefer my indexes to start at:

    Results (126 votes). Check out past polls.