in reply to RFC: 100 PDL Exercises (ported from numpy)

PM is a good choice for this project. It's got the discussion of the best solution(s) to the question. Hopefully someone jumps in an offers a better solution than mine. It feels a little ugly explicitly stating the indices.

Given two arrays, X and Y, construct the Cauchy matrix C (Cij =1/(xi - yj))

First, what's a Cauchy matrix? Ahh, this question is just how to construct a matrix from 2 arrays, where no elements from one array are in the other. Just create a sequence of number for the first array and then make the second array 0.5 more than the first array to get the inputs. Here's a brute force method. 
use PDL;
use PDL::NiceSlice;

my $x = sequence(8);
my $y = $x + 0.5;

my ($nx, $ny) = (nelem($x), nelem($y));

my $C = zeroes($nx, $ny);

for (my $i = 0; $i < $nx; $i++) {
    for (my $j = 0; $j < $ny; $j++) {
        $C($i,$j) .= 1/($x($i) - $y($j));
    }
}

print $C;

[download]
I like the PDL::NiceSlice for indexing. It makes sense to me. I could have also created the matrix with  my $C = outer($x, $y); or gotten the size of the arrays with  $x->getdim(0) and if I grokked threading rather than just skimming PDL::Threading, this might look way cooler.

NB the ".=" in the assignment breaks the link between the matrix and the 2 arrays. It's important.

That was just the Cauchy matrix. Some people want the Cauchy determinant (as long as the 2 arrays are the same size). Easy! Just import PDL::MatrixOps 

use PDL::MatrixOps;

print det $C;

[download]

Sometimes I can think of 6 impossible LDAP attributes before breakfast.

YAPC::Europe::2018 — Hmmm, need to talk to work about sending me.

Replies are listed 'Best First'.
Re: 47. Construct the Cauchy matrix
by Ea (Chaplain) on Feb 20, 2020 at 18:03 UTC

    This is Much Better

    When I was young, I coded in for loops, but when I read the documentation, I put away these childish things and learned to Thread Broadcast 
    $x = sequence(8);
$y = sequence(7) + 0.5;       
 
$c = 1/($x->dummy(0,$y->nelem)->transpose - $y->dummy(0,$x->nelem));

    [download]
    TaDAAA!!

    You create 2 vectors (of different sizes), inflate it along the other dimension using dummy to fit the other vector's size (using nelem), flip one of them using transpose and do the calculation in one line. PDL takes care of the loops and does it faster than you can in Perl.

    I had to transpose $x to get the same result as the for loops above, in order to prove they are the same by pdl> p $C - $c

    Ea

    Sometimes I can think of 6 impossible LDAP attributes before breakfast.

    Mojoconf was great!

[reply]
[d/l]
[select]
      Excellent use of broadcasting (the feature formerly and confusingly known as "threading")! To make that broadcast even better, you'd dummy not by the nelem, but explicitly with dim(1) (or 0) as someone else on here did.
[reply]
[d/l]
[select]

In Section Meditations