A PyTorch implementation of Algikar et.al.'s projection statistics
Last, we provide the code of the main parts of our translation of the Matlab code of Algikar et.al., 2023 in the following, for you to be able to check its correctness.
```
def projection_statistics(H: Tensor) -> Tensor:
"""
Args:
H: (n x d)-dim Tensor, i.e. number of data points x dimensionality. NOTE:
in the original code, this is taken to be X with an appended ones column.
Returns:
A (n)-dim Tensor of projection statistics.
"""
dtype = H.dtype
device = H.device
m, n = H.shape
M = torch.median(H, dim=0, keepdim=True).values # (1 x d) row vector
u = torch.zeros(m, n, dtype=dtype, device=device) # i.e. (n x d) matrix
v = torch.zeros(m, n, dtype=dtype, device=device)
z = torch.zeros(m, 1, dtype=dtype, device=device) # i.e. (n x 1) matrix
P = torch.zeros(m, m, dtype=dtype, device=device) # i.e. (n x n) matrix
eps = 1e-6 # avoiding divide-by-zero issues
for kk in range(m):
u[kk, :] = H[kk, :] - M # looping over data points
v[kk, :] = u[kk, :] / max(torch.linalg.norm(u[kk, :]), eps)
for ii in range(m):
z[ii, :] = torch.dot(H[ii, :], v[kk, :])
zmed = torch.median(z, dim=0, keepdim=True).values
MAD = 1.4826 * (1 + (15 / (m))) * torch.median(torch.abs(z - zmed))
for ii in range(m):
P[kk, ii] = torch.abs(z[ii] - zmed) / max(MAD, eps)
PS = torch.amax(P, dim=0)
return PS
def _compute_projection_statistics_weights(H: Tensor, PS: Tensor) -> Tensor:
"""This computes the weights for the projection statistics, according to the procedure
in the original Matlab code:
https://github.com/apooja1/GP-Huber/blob/fee038963b471eb198d59b22d57b89e69f451d8c/Experiments/Friedman
"""
niu = torch.sum(H != 0, dim=-1)
cutoff_PS = torch.zeros_like(niu)
for i in range(len(niu)):
# for chi2, see this post for correspondence of ppf with invcdf:
# https://stackoverflow.com/questions/53019080/chi2inv-in-python
# the 0.975 was copied from the original matlab code:
# https://github.com/apooja1/GP-Huber/blob/fee038963b471eb198d59b22d57b89e69f451d8c/Experiments/Friedman.m#L186
cutoff_PS[i] = chi2.ppf(0.975, df=niu[i].item())
weights = (cutoff_PS / PS.square()).clamp(max=1.0)
return weights
def compute_projection_statistics_weights(X: Tensor) -> Tensor:
# X: n x d
n, d = X.shape
H = torch.cat((X, torch.ones(n, 1)), dim=-1) # n x (d + 1)
PS = projection_statistics(H) # (n x d) -> n
return _compute_projection_statistics_weights(H, PS)
```
For the Huber loss, we use $\epsilon = 0.45$ and $b = 0.5$, identical to the setting of the [Github repository](https://github.com/apooja1/GP-Huber/blob/fee038963b471eb198d59b22d57b89e69f451d8c/basics/lik_huber.m#L191) of Algikar et.al., 2023. All reported results were generated with 32 independent replications, reporting the mean and standard error.
An additional implementation detail: the `median` and `max` calls take the extra `dim=0` argument in PyTorch, as the Matlab function are defined to compute the row-wise maximum. Also, we added a small numerical constant `eps` in the code of `projection_statistics` to avoid division-by-zero errors. The rest of the code should be an almost verbatim translation.
**A word of thanks** Thank you again for your valuable suggestions, we believe they make the paper stronger. Let us know if you have additional questions or if all your remaining concerns have been addressed.