Dimitrios@0: function [W,H,errs,vout] = nmf_beta(V,r,varargin) Dimitrios@0: % function [W,H,errs,vout] = nmf_beta(V,r,varargin) Dimitrios@0: % Dimitrios@0: % Implements NMF using the beta-divergence [1]: Dimitrios@0: % Dimitrios@0: % min D(V||W*H) s.t. W>=0, H>=0 Dimitrios@0: % Dimitrios@0: % / Dimitrios@0: % | sum(V(:).^beta + (beta-1)*R(:).^beta - ... Dimitrios@0: % | beta*V(:).*R(:).^(beta-1)) / ... Dimitrios@0: % | (beta*(beta-1)) (beta \in{0 1} Dimitrios@0: % where D(V||R) = | Dimitrios@0: % | sum(V(:).*log(V(:)./R(:)) - V(:) + R(:)) (beta=1) Dimitrios@0: % | Dimitrios@0: % | sum(V(:)./R(:) - log(V(:)./R(:)) - 1) (beta=0) Dimitrios@0: % \ Dimitrios@0: % Dimitrios@0: % This divergence reduces to the following interesting distances for Dimitrios@0: % certain values of beta: Dimitrios@0: % Dimitrios@0: % - Itakura-Saito (beta=0) Dimitrios@0: % - I-divergence (beta=1) Dimitrios@0: % - Euclidean distance (beta=2) Dimitrios@0: % Dimitrios@0: % Inputs: (all except V and r are optional and passed in in name-value pairs) Dimitrios@0: % V [mat] - Input matrix (n x m) Dimitrios@0: % r [num] - Rank of the decomposition Dimitrios@0: % beta [num] - beta parameter [0] Dimitrios@0: % niter [num] - Max number of iterations to use [100] Dimitrios@0: % thresh [num] - Number between 0 and 1 used to determine convergence; Dimitrios@0: % the algorithm has considered to have converged when: Dimitrios@0: % (err(t-1)-err(t))/(err(1)-err(t)) < thresh Dimitrios@0: % ignored if thesh is empty [[]] Dimitrios@0: % norm_w [num] - Type of normalization to use for columns of W [1] Dimitrios@0: % can be 0 (none), 1 (1-norm), or 2 (2-norm) Dimitrios@0: % norm_h [num] - Type of normalization to use for rows of H [0] Dimitrios@0: % can be 0 (none), 1 (1-norm), 2 (2-norm), or 'a' (sum(H(:))=1) Dimitrios@0: % verb [num] - Verbosity level (0-3, 0 means silent) [1] Dimitrios@0: % W0 [mat] - Initial W values (n x r) [[]] Dimitrios@0: % empty means initialize randomly Dimitrios@0: % H0 [mat] - Initial H values (r x m) [[]] Dimitrios@0: % empty means initialize randomly Dimitrios@0: % W [mat] - Fixed value of W (n x r) [[]] Dimitrios@0: % empty means we should update W at each iteration while Dimitrios@0: % passing in a matrix means that W will be fixed Dimitrios@0: % H [mat] - Fixed value of H (r x m) [[]] Dimitrios@0: % empty means we should update H at each iteration while Dimitrios@0: % passing in a matrix means that H will be fixed Dimitrios@0: % myeps [num] - Small value to add to denominator of updates [1e-20] Dimitrios@0: % Dimitrios@0: % Outputs: Dimitrios@0: % W [mat] - Basis matrix (n x r) Dimitrios@0: % H [mat] - Weight matrix (r x m) Dimitrios@0: % errs [vec] - Error of each iteration of the algorithm Dimitrios@0: % Dimitrios@0: % [1] Cichocki, A. and Amari, S.I. and Zdunek, R. and Kompass, R. and Dimitrios@0: % Hori, G. and He, Z. Extended SMART Algorithms for Non-negative Matrix Dimitrios@0: % Factorization, Artificial Intelligence and Soft Computing, 2006 Dimitrios@0: % Dimitrios@0: % 2010-01-14 Graham Grindlay (grindlay@ee.columbia.edu) Dimitrios@0: Dimitrios@0: % Copyright (C) 2008-2028 Graham Grindlay (grindlay@ee.columbia.edu) Dimitrios@0: % Dimitrios@0: % This program is free software: you can redistribute it and/or modify Dimitrios@0: % it under the terms of the GNU General Public License as published by Dimitrios@0: % the Free Software Foundation, either version 3 of the License, or Dimitrios@0: % (at your option) any later version. Dimitrios@0: % Dimitrios@0: % This program is distributed in the hope that it will be useful, Dimitrios@0: % but WITHOUT ANY WARRANTY; without even the implied warranty of Dimitrios@0: % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the Dimitrios@0: % GNU General Public License for more details. Dimitrios@0: % Dimitrios@0: % You should have received a copy of the GNU General Public License Dimitrios@0: % along with this program. If not, see . Dimitrios@0: Dimitrios@0: % do some sanity checks Dimitrios@0: if min(min(V)) < 0 Dimitrios@0: error('Matrix entries can not be negative'); Dimitrios@0: end Dimitrios@0: if min(sum(V,2)) == 0 Dimitrios@0: error('Not all entries in a row can be zero'); Dimitrios@0: end Dimitrios@0: Dimitrios@0: [n,m] = size(V); Dimitrios@0: Dimitrios@0: % process arguments Dimitrios@0: [beta, niter, thresh, norm_w, norm_h, verb, myeps, W0, H0, W, H] = ... Dimitrios@0: parse_opt(varargin, 'beta', 0, 'niter', 100, 'thresh', [], ... Dimitrios@0: 'norm_w', 1, 'norm_h', 0, 'verb', 1, 'myeps', 1e-20, ... Dimitrios@0: 'W0', [], 'H0', [], 'W', [], 'H', []); Dimitrios@0: Dimitrios@0: % initialize W based on what we got passed Dimitrios@0: if isempty(W) Dimitrios@0: if isempty(W0) Dimitrios@0: W = rand(n,r); Dimitrios@0: else Dimitrios@0: W = W0; Dimitrios@0: end Dimitrios@0: update_W = true; Dimitrios@0: else Dimitrios@0: update_W = false; Dimitrios@0: end Dimitrios@0: Dimitrios@0: % initialize H based on what we got passed Dimitrios@0: if isempty(H) Dimitrios@0: if isempty(H0) Dimitrios@0: H = rand(r,m); Dimitrios@0: else Dimitrios@0: H = H0; Dimitrios@0: end Dimitrios@0: update_H = true; Dimitrios@0: else % we aren't H Dimitrios@0: update_H = false; Dimitrios@0: end Dimitrios@0: Dimitrios@0: if norm_w ~= 0 Dimitrios@0: % normalize W Dimitrios@0: W = normalize_W(W,norm_w); Dimitrios@0: end Dimitrios@0: Dimitrios@0: if norm_h ~= 0 Dimitrios@0: % normalize H Dimitrios@0: H = normalize_H(H,norm_h); Dimitrios@0: end Dimitrios@0: Dimitrios@0: % initial reconstruction Dimitrios@0: R = W*H; Dimitrios@0: Dimitrios@0: errs = zeros(niter,1); Dimitrios@0: for t = 1:niter Dimitrios@0: % update W if requested Dimitrios@0: if update_W Dimitrios@0: W = W .* ( ((R.^(beta-2) .* V)*H') ./ max(R.^(beta-1)*H', myeps) ); Dimitrios@0: if norm_w ~= 0 Dimitrios@0: W = normalize_W(W,norm_w); Dimitrios@0: end Dimitrios@0: end Dimitrios@0: Dimitrios@0: % update reconstruction Dimitrios@0: R = W*H; Dimitrios@0: Dimitrios@0: % update H if requested Dimitrios@0: if update_H Dimitrios@0: H = H .* ( (W'*(R.^(beta-2) .* V)) ./ max(W'*R.^(beta-1), myeps) ); Dimitrios@0: if norm_h ~= 0 Dimitrios@0: H = normalize_H(H,norm_h); Dimitrios@0: end Dimitrios@0: end Dimitrios@0: Dimitrios@0: % update reconstruction Dimitrios@0: R = W*H; Dimitrios@0: Dimitrios@0: % compute beta-divergence Dimitrios@0: switch beta Dimitrios@0: case 0 Dimitrios@0: errs(t) = sum(V(:)./R(:) - log(V(:)./R(:)) - 1); Dimitrios@0: case 1 Dimitrios@0: errs(t) = sum(V(:).*log(V(:)./R(:)) - V(:) + R(:)); Dimitrios@0: case 2 Dimitrios@0: errs(t) = sum(sum((V-W*H).^2)); Dimitrios@0: otherwise Dimitrios@0: errs(t) = sum(V(:).^beta + (beta-1)*R(:).^beta - beta*V(:).*R(:).^(beta-1)) / ... Dimitrios@0: (beta*(beta-1)); Dimitrios@0: end Dimitrios@0: Dimitrios@0: % display error if asked Dimitrios@0: if verb >= 3 Dimitrios@0: disp(['nmf_beta: iter=' num2str(t) ', err=' num2str(errs(t)) ... Dimitrios@0: '(beta=' num2str(beta) ')']); Dimitrios@0: end Dimitrios@0: Dimitrios@0: % check for convergence if asked Dimitrios@0: if ~isempty(thresh) Dimitrios@0: if t > 2 Dimitrios@0: if (errs(t-1)-errs(t))/(errs(1)-errs(t-1)) < thresh Dimitrios@0: break; Dimitrios@0: end Dimitrios@0: end Dimitrios@0: end Dimitrios@0: end Dimitrios@0: Dimitrios@0: % display error if asked Dimitrios@0: if verb >= 2 Dimitrios@0: disp(['nmf_beta: final_err=' num2str(errs(t)) '(beta=' num2str(beta) ')']); Dimitrios@0: end Dimitrios@0: Dimitrios@0: % if we broke early, get rid of extra 0s in the errs vector Dimitrios@0: errs = errs(1:t); Dimitrios@0: Dimitrios@0: % needed to conform to function signature required by nmf_alg Dimitrios@0: vout = {};