Consistent Data Matching via Laplacian Clusters

# Illustrative sketch (not from the paper)
import torch
import torch.nn.functional as F
from sklearn.cluster import KMeans, SpectralClustering

# X, Y: point clouds (N x D) tensors
X = torch.randn(100, 3)
Y = torch.randn(120, 3)

# 1. similarity graphs via RBF kernel
pairwise_xy = torch.cdist(X, X, p=2)
Sx = torch.exp(-pairwise_xy**2 / (2 * 1.0**2))  # sigma=1.0
pairwise_yy = torch.cdist(Y, Y, p=2)
Sy = torch.exp(-pairwise_yy**2 / (2 * 1.0**2))

# 2. Laplacian matrices
Lx = torch.diag(Sx.sum(dim=1)) - Sx
Ly = torch.diag(Sy.sum(dim=1)) - Sy

# 3. Laplacian OT (Sinkhorn with Laplacian regularization placeholder)
#   min_{T} <C, T> + lambda/2 * (Tr(T^T Lx T) + Tr(T Ly T^T))
#   subject to T 1 = a, T^T 1 = b, T >= 0
C = torch.cdist(X, Y, p=2)  # cost matrix
lambda_reg = 0.1
# Simple Sinkhorn iteration (without Laplacian term for brevity)
a = torch.full((X.shape[0],), 1.0 / X.shape[0])
b = torch.full((Y.shape[0],), 1.0 / Y.shape[0])
K = torch.exp(-C / lambda_reg)
U = K.clone()
for _ in range(100):
    U = U * (a / (U @ b)).unsqueeze(1)
    U = U * (b / (U.t() @ a)).unsqueeze(0)
T = U  # coupling matrix

# 4. Refined Simultaneous Clustering (RSC) on coupling matrix
k = 5  # number of clusters (example)
km = KMeans(n_clusters=k).fit(T.detach().cpu().numpy())
switch_X = torch.tensor(km.labels_).unsqueeze(1) == torch.arange(k)
switch_Y = torch.tensor(km.labels_).unsqueeze(0) == torch.arange(k)
# switch matrices are binary masks indicating cluster membership

# 5. Refine similarity matrices
Sx_refined = Sx * switch_X.float()
Sy_refined = Sy * switch_Y.float()

# 6. Updated Laplacians and spectral clustering
Lx_ref = torch.diag(Sx_refined.sum(dim=1)) - Sx_refined
Ly_ref = torch.diag(Sy_refined.sum(dim=1)) - Sy_refined
sc = SpectralClustering(n_clusters=k, affinity='precomputed')
clusters_X = sc.fit_predict(Sx_refined.detach().cpu().numpy())
clusters_Y = sc.fit_predict(Sy_refined.detach().cpu().numpy())

// Illustrative sketch (not from the paper)
const tf = require('@tensorflow/tfjs-node');
const { KMeans } = require('ml-kmeans'); // placeholder library
const { SpectralClustering } = require('ml-spectral-clustering'); // placeholder

// X, Y: point clouds as tensors (N x D)
const X = tf.randomNormal([100, 3]);
const Y = tf.randomNormal([120, 3]);

// 1. similarity graphs via RBF kernel
function rbfSimilarity(A) {
  const dists = tf.norm(tf.expandDims(A, 1).sub(tf.expandDims(A, 0)), 'euclidean', 2);
  return tf.exp(tf.neg(tf.square(dists).div(2 * Math.pow(1.0, 2)))); // sigma=1.0
}
const Sx = rbfSimilarity(X);
const Sy = rbfSimilarity(Y);

// 2. Laplacian matrices
function laplacian(S) {
  const degree = tf.sum(S, 1);
  return tf.diag(degree).sub(S);
}
const Lx = laplacian(Sx);
const Ly = laplacian(Sy);

// 3. Laplacian OT (Sinkhorn placeholder)
const C = tf.norm(tf.expandDims(X, 1).sub(tf.expandDims(Y, 0)), 'euclidean', 2);
const lambdaReg = 0.1;
const a = tf.fill([X.shape[0]], 1.0 / X.shape[0]);
const b = tf.fill([Y.shape[0]], 1.0 / Y.shape[0]);
let K = tf.exp(tf.neg(C.div(lambdaReg)));
let U = K.clone();
for (let i = 0; i < 100; i++) {
  const colSum = tf.matMul(U, b.expandDims(1)).squeeze();
  U = U.mul(a.div(colSum).expandDims(1));
  const rowSum = tf.matMul(U.transpose(), a.expandDims(1)).squeeze();
  U = U.mul(b.div(rowSum).expandDims(0));
}
const T = U; // coupling matrix

// 4. Refined Simultaneous Clustering (RSC) on coupling matrix
const k = 5; // number of clusters (example)
const km = new KMeans(T.arraySync(), { k });
const labels = km.clusters; // array of length N+M
const switchX = tf.tensor(labels.slice(0, X.shape[0])).expandDims(1).equal(tf.range(k));
const switchY = tf.tensor(labels.slice(X.shape[0])).expandDims(0).equal(tf.range(k));

// 5. Refine similarity matrices
const SxRefined = Sx.mul(switchX.cast('float32'));
const SyRefined = Sy.mul(switchY.cast('float32'));

// 6. Updated Laplacians and spectral clustering
const LxRef = laplacian(SxRefined);
const LyRef = laplacian(SyRefined);
const scX = new SpectralClustering({ nClusters: k, affinity: 'precomputed' });
const clustersX = scX.fitPredict(SxRefined.arraySync());
const scY = new SpectralClustering({ nClusters: k, affinity: 'precomputed' });
const clustersY = scY.fitPredict(SyRefined.arraySync());