47. Construct the Cauchy matrix

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 (C_ij =1/(x_i - y_j))

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.

Comment on 47. Construct the Cauchy matrix Select or Download Code

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]
Re^2: 47. Construct the Cauchy matrix by etj (Deacon) on Mar 16, 2022 at 21:05 UTC
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]


P is for Practical
	PerlMonks

47. Construct the Cauchy matrix

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

This is Much Better

Ea

Given two arrays, X and Y, construct the Cauchy matrix C (C_ij =1/(x_i - y_j))