File size: 4,740 Bytes
80a72c3 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 | import numpy as np
def mse_loss_at_residue(V_hat, Y, residue_index):
# loss contribution of residue given current assignment
# assumes Y is symmetric
d = residue_index
# avoid double-counting the diagonal
return 2*((V_hat[d]@V_hat.T - Y[d])**2).sum() - (V_hat[d]@V_hat[d] - Y[d,d])**2
def null_mse_loss_at_residue(Y, residue_index):
# loss contribution of residue when unassigned
d = residue_index
return 2*(Y[d]**2).sum() - Y[d,d]**2
def mse_loss(V_hat, Y):
# only use to compute initial loss
return ((V_hat@V_hat.T - Y)**2).sum()
def nll_loss(V_hat, Y):
"""We now interpret ys as probabilities.
We want to minimise the negative log probability.
NB Y must be symmetric and [0-1]
"""
# only use to compute initial loss
log_pij = (V_hat@V_hat.T) * np.log(Y) + (1-(V_hat@V_hat.T))*np.log(1-Y)
# simple sum double counts off diagonals but single counts on diagonals, so need
# to add diagonal again before dividing by 2
return -0.5*(log_pij.sum() + np.diag(log_pij).sum())
def nll_loss_at_residue(V_hat, Y, residue_index):
# d is residue_index
d = residue_index # V_hat[d]@V_hat.T is just Adj without instantiating A
# n.b. whereas in case of mse we double-count off-diagonals, here we don't, hence the difference
log_pdj = (V_hat[d]@V_hat.T) * np.log(Y[d]) + (1-(V_hat[d]@V_hat.T))*np.log(1-Y[d])
return -log_pdj.sum()
def null_nll_loss_at_residue(Y, residue_index):
return - np.log(1-Y[residue_index]).sum()
def greedy_V(Y, N_iters=3, K_init=4, cost_type="mse"):
"""
Learn a binary matrix V, with at most one nonzero entry per row, to minimize
|| VV' - Y ||_2
This is done by initializing a V of all zeros, and then doing a greedy optimization.
V is initially D x K_init, where D is the number of residues (i.e., where Y is DxD).
K is learned automatically; different values of K_init will not change the result but
might be more or less efficient in terms of memory usage. The code keeps track of the
current number of nonzero columns of V, and when all the columns are full it adds extra
columns of zeros at the end.
Each iteration sweeps through all residues (i.e. rows of V) once.
The implementation relies on additivity of loss to reduce computational cost
by only computing required increments to the loss at each iteration.
INPUTS:
Y: matrix of predictions with entries in [0, 1]
N_iters: number of iterations
K_init: initial number of columns of V, adjust this for tweaking performance
"""
Y = (Y + Y.T) / 2 # required for consistency with model_v1
V_hat = np.zeros((Y.shape[0], K_init), dtype=np.uint8)
if cost_type == "mse":
loss = mse_loss(V_hat, Y) # initial loss for the zero matrix V_hat = 0
elif cost_type == "nll":
loss = nll_loss(V_hat, Y)
else:
raise ValueError(cost_type)
K_max = K_init # track number of columns in K_max
for it in range(N_iters):
for d in range(V_hat.shape[0]):
# d is a residue index
# compute the loss, excluding contribution of this residue
# TODO check minus signs are consistent
loss_minus_d = loss - mse_loss_at_residue(V_hat, Y, d)
# sweep through all K+1 options and compute what the contribution to the loss would be
V_hat[d] *= 0
# loss with no assignment
L0 = loss_minus_d + null_mse_loss_at_residue(Y, d)
L_opt = np.zeros(K_max)
for k in range(K_max):
# note this could be vectorized with a bit of work
V_hat[d,k] = 1
L_opt[k] = loss_minus_d + mse_loss_at_residue(V_hat, Y, d)
V_hat[d,k] = 0
# select the option which minimizes the squared error
z = np.argmin(L_opt)
if L_opt[z] < L0:
V_hat[d,z] = 1
# update loss
loss = loss_minus_d + mse_loss_at_residue(V_hat, Y, d)
if z == K_max-1:
# Expand V_hat to make room for extra potential clusters (no nonzero columns remain)
# TODO: note that if a cluster is "removed" on a later iteration we will miss this.
# probably doesn't matter, just means it is occasionally possible that the V
# that is returned could sometimes have a few nonzero columns.
V_hat = np.concatenate((V_hat, np.zeros_like(V_hat)), -1)
K_max = V_hat.shape[1]
# drop columns that are zeros (i.e. unused clusters) before returning
empty = V_hat.sum(0) == 0
return V_hat[:,~empty], loss
|