Distance matrix computation

[bcsj]

I'm interesting in computing the distance matrix for a very large dataset. Now, there is a formula,

XXt = X * X^T

dist[i,j]^2 = Z1 - 2*XXt + Z2

where Z1 adds to each row the vector of the diagonal elements of XXt, and Z2 adds adds to each column the vector of the diagonal elements of XXt.

Is there a good way in slate to extract the diagonal of a matrix, and following that to broadcast-add a vector to each row/column of a matrix?

I've been looking through the example codes in github and trying to search the documentation on bitbucket, but I don't feel like I've gotten closer to a good way to do either so far.

[Dalal Sukkari]

What is your matrix type? i.e, Is it Matrix, TriangularBandMatrix, HermitianBandMatrix, ...?

In src/heev.cc we copy the diagonal and super-diagonal of a distributed matrix by calling:

// Copy diagonal and super-diagonal to vectors.                                                              
internal::copyhb2st(Aband, Lambda, E);

where Lambda will have the diagonal of the matrix Aband and E will have the super-diagonal.

Similarly, we call internal::copytb2bd(Aband, Sigma, E) in src/gesvd.cc.

Hopefully that helps you to do your computations.

--
You received this message because you are subscribed to the Google Groups "SLATE User" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slate-user+...@icl.utk.edu.
To view this discussion on the web visit https://groups.google.com/a/icl.utk.edu/d/msgid/slate-user/e3199863-5d0a-4b3f-811a-23af9c0010den%40icl.utk.edu.

[bcsj]

The stuff in src/heev.cc may indeed be something like what I'm looking for. I'll have to dig a bit more around in the source code I think.

Not exactly sure what you are asking for about the matrix type, but I can try to elaborate:

I have a dataset represented by a matrix X, slate::Matrix<scalar_t> X(...); where scalar_t is say float or double.

There is no inherent structure to X, it is dense and very tall, each row being a data sample with a number of features.

I'm interested in computing the distance matrix D which has elements d[i,j] = ||x[i,:] - x[j,:]||. 

This could be done in various ways:

 - 1. I could loop through the data a bunch

 - 2. I could divide D into different tiles and distribute the computation on each tile

 - 3. I could do the above, but write a GPU kernel for taking care of the actual computations

 - 4. I could use a formula involving a matrix product, 2 row/columns-wise additions and a element-wise squareroot

There are possibly more ways. 

1. is not very efficient, 

2. and 3. are promising, likely the best given the required skillset, but will require a lot of code optimization to get good performance, 

4. lets me use the existing gemm implementation and hopefully skip a lot of optimizing myself.

So I'm hoping to do 4. I might attempt 2.&3. on the side, but I'm also asked to try to get something working pretty fast, so I can't just take my time on 2.&3.

The formula for 4. requires computing the matrix product X * X^T, that matrix is of course symmetric, and it is the diagonal of that resulting matrix I'd want to extract.

[bcsj]

I've been poking around in the code in heev.cc and the code for he2hbGather(...) in HermitianBandMatrix.hh and I think I've managed to cook up some code of my own that does a similar thing.

In doing this, I've noticed that sometimes I may have the tiles of a matrix available on my mpi_rank, despite the tile being owned by another mpi_rank. In other words, if I check A.tileIsLocal(i, j) and get False, also when a local version of the tile is available. Is there a way to check if a tile is available? i.e. something like A.tileIsLocallyAvailable(i, j)?

[bcsj]

Okay, I think I managed to find something which checks if a non-local tile is available locally: A.tileExists(i, j) seems to do that, assuming I understand it correctly.

[Mark Gates]

Yes, tileExists( i, j, device=Host ) tells if a tile is present, whether it is local or remote.

If you're using GPUs, it's possible for a tile to exist but not be valid, e.g., the tile exists on both host and device, and the device copy is Modified, so the host copy is Invalid. It's best to do a tileGetForReading( i, j, device=Host ) beforehand to ensure that the host copy is valid. (This will hopefully be simplified in future versions.)

Mark

[bcsj]

Great, I did notice the tileGetForReading being used in various places, but I didn't understand what it was doing.

Is there a method for getting the global index (p,q) of a matrix entry from the tile index (ii, jj) and tile-local index (i,j)? 

Or do I need to compute it myself?

[Mark Gates]

On Thu, Jun 15, 2023 at 4:38 AM bcsj <bjornje...@gmail.com> wrote:

Is there a method for getting the global index (p,q) of a matrix entry from the tile index (ii, jj) and tile-local index (i,j)? 

Or do I need to compute it myself?

I don't there is. If you assume fixed block size nb, it should be easy:

    p = ii*nb + i

(although I would usually use i for tile index, ii for index within tile, and p, q are usually the MPI process grid dimension).

Related, if you assume 2D block-cyclic ScaLAPACK layout, there are local2global and global2local. They are marked as internal routines because we often don't use ScaLAPACK layout.

Mark

--

Innovative Computing Laboratory

University of Tennessee, Knoxville