Beefy Boxes and Bandwidth Generously Provided by pair Networks
Clear questions and runnable code
get the best and fastest answer

47. Construct the Cauchy matrix

by Ea (Hermit)
on May 18, 2018 at 15:32 UTC ( #1214858=note: print w/replies, xml ) Need Help??

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;
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;

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 (Hermit) 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
    $x = sequence(8); $y = sequence(7) + 0.5; $c = 1/($x->dummy(0,$y->nelem)->transpose - $y->dummy(0,$x->nelem));

    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


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

    Mojoconf was great!

Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://1214858]
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others pondering the Monastery: (4)
As of 2020-03-29 10:08 GMT
Find Nodes?
    Voting Booth?
    To "Disagree to disagree" means to:

    Results (169 votes). Check out past polls.