Commit b97d94c6 authored by Maxime's avatar Maxime

Code for beta_div != 1

parent dd9fedb2
......@@ -4,7 +4,7 @@ Sparse NMF
Example script for Sparse NMF with beta-divergence distortion function and L1 penalty on the activations.
# Disclaimer
This port has not been fully tested, and so far only the algorithm only runs with `beta_div=1`.
This port has not been fully tested, especially with parameters `beta_div != 1`.
# Copyright notice and License
Copyright (C) 2015 Mitsubishi Electric Research Labs
......
......@@ -96,7 +96,7 @@ def sparse_nmf(v, beta=1, cf='kl', sparsity=0, max_iter=100, conv_eps=0,
if init_h is None: ## init h
h = np.random.rand(r,n)
elif init_h=='ones': # what the hell is that??
elif init_h=='ones':
print('INFO: initializing H with ones')
h = np.ones((r,n))
else:
......@@ -111,20 +111,16 @@ def sparse_nmf(v, beta=1, cf='kl', sparsity=0, max_iter=100, conv_eps=0,
if type(sparsity) in (int, float): # sparsity per matrix entry
sparsity = np.ones((r,n))*sparsity
elif sparsity.shape[1] == 1:
sparsity = repmat(sparsity, 1, n) ## to debug, use numpy
sparsity = np.tile(sparsity, (1,n)) ## to debug
## Normalize the columns of W and rescale H accordingly
wn = ((w**2).sum(axis=0))**0.5
w = w/np.tile(wn, (w.shape[0], 1)) ## surely wrong
h = h*np.tile(wn, (h.shape[1], 1)).T
##wn = sqrt(sum(w.^2));
##w = bsxfun(@rdivide,w,wn);
##h = bsxfun(@times, h,wn');
wn = ((w**2).sum(axis=0))**0.5 ##wn = sqrt(sum(w.^2));
w = w/np.tile(wn, (w.shape[0], 1)) ##w = bsxfun(@rdivide,w,wn);
h = h*np.tile(wn, (h.shape[1], 1)).T ##h = bsxfun(@times, h,wn');
## Internal parameters
flr = 1e-9
lamb = w.dot(h) # used to be lambda
lamb = w.dot(h) # used to be lambda (but lambda is a reserved keyword in Python)
lamb[lamb<flr]=flr
last_cost = np.inf
......@@ -153,16 +149,16 @@ def sparse_nmf(v, beta=1, cf='kl', sparsity=0, max_iter=100, conv_eps=0,
#dmh = w(:, h_ind)' * (v ./ lamb);
#h(h_ind, :) = bsxfun(@rdivide, h(h_ind, :) .* dmh, dph);
elif div_beta==2:
dph = w[:, h_ind].T * lamb +sparsity
dph = np.dot(w[:, h_ind].T, lamb) +sparsity
dph[dph<flr]=flr
dmh = np.multiply(w[:, h_ind].T, v)
dmh = np.dot(w[:, h_ind].T, v)
#dph = w(:, h_ind)' * lambda + params.sparsity;
#dph = max(dph, flr);
#dmh = w(:, h_ind)' * v;
else:
dph = w[:, h_ind].T * lamb**(div_beta-1) +sparsity
dph = np.dot(w[:, h_ind].T, lamb**(div_beta-1)) +sparsity
dph[dph<flr]=flr
dmh = np.multiply(w[:, h_ind].T, v*lamb**(div_beta-2))
dmh = np.dot(w[:, h_ind].T, v*lamb**(div_beta-2))
h[h_ind, :] = h[h_ind,:]*dmh/dph
#dph = w(:, h_ind)' * lambda.^(div_beta - 1) + params.sparsity;
#dph = max(dph, flr);
......@@ -191,20 +187,38 @@ def sparse_nmf(v, beta=1, cf='kl', sparsity=0, max_iter=100, conv_eps=0,
w[:, w_ind] = w[:, w_ind] * dmw /dpw
#w(:, w_ind) = w(:,w_ind) .* dmw ./ dpw;
elif div_beta==2:
pass
dpw1 = np.dot(lamb , h[w_ind, :].T)
dpw2 = (np.dot(v, h[w_ind,:].T)*w[:,w_ind]).sum(0)
dpw2 = np.tile(dpw2, (w.shape[0], 1)) * w[:, w_ind]
dpw = dpw1 + dpw2
dpw[dpw<flr] = flr
#dpw = lamb * h(w_ind, :)' ...
# + bsxfun(@times, sum(v * h(w_ind, :)' .* w(:, w_ind)), w(:, w_ind));
#dpw = max(dpw, flr);
dmw1 = np.dot(v, h[w_ind,:].T)
dmw2 = (np.dot(lamb, h[w_ind,:].T)*w[:,w_ind]).sum(0)
dmw2 = np.tile(dmw2, (w.shape[0], 1)) * w[:, w_ind]
dmw = dmw1 + dmw2
#dmw = v * h(w_ind, :)' + ...
# bsxfun(@times, sum(lambda * h(w_ind, :)' .* w(:, w_ind)), w(:, w_ind));
w[:,w_ind] = w[:,w_ind] * dmw / dpw
#w(:, w_ind) = w(:,w_ind) .* dmw ./ dpw;
else:
pass
dpw1 = np.dot(lamb**(div_beta-1), h[w_ind, :].T)
dpw2 = (np.dot(v*lamb**(div_beta-2), h[w_ind,:].T) * w[:, w_ind]).sum(0)
dpw2 = np.tile(dpw2, (w.shape[0], 1)) * w[:, w_ind]
dpw = dpw1 + dpw2
dpw[dpw<flr] = flr
#dpw = lamb.^(div_beta - 1) * h(w_ind, :)' ...
# + bsxfun(@times, ...
# sum((v .* lambda.^(div_beta - 2)) * h(w_ind, :)' .* w(:, w_ind)), ...
# w(:, w_ind));
#dpw = max(dpw, flr);
dmw1 = np.dot(v*lamb**(div_beta-2), h[w_ind, :].T)
dmw2 = (np.dot(lamb**(div_beta-1), h[w_ind, :].T) * w[:, w_ind]).sum(0)
dmw2 = np.tile(dmw2, (w.shape[0], 1)) * w[:, w_ind]
dmw = dmw1 + dmw2
w[:,w_ind] = w[:,w_ind] * dmw / dpw
#dmw = (v .* lamb.^(div_beta - 2)) * h(w_ind, :)' ...
# + bsxfun(@times, ...
# sum(lambda.^(div_beta - 1) * h(w_ind, :)' .* w(:, w_ind)), w(:, w_ind));
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment