MVN optimization with Cholesky + minor additions

README.md

@@ -31,7 +31,8 @@ Contact: [Reshad Hosseini](mailto:reshad.hosseini@ut.ac.ir) or [Mohamadreza Mash
    |--- manopt/           - manifold optimization toolbox by N. Boumal and B. Mishra
    |--- matlab-xunit/     - xUnit tests for Matlab
 
- </pre>
+</pre>
+
 ## HOW TO CITE ##
 
 If you are using this toolbox in your research please cite the following paper:
@@ -89,12 +90,14 @@ Using MixEst, you can fit arbitrary distributions (usually mixture models) to yo
 <section>
 <img src="http://visionlab.ut.ac.ir/mixest/docs/examples/img/example3.gif">
 </section>
+
 An example showing the toolbox fitting a mixture on von Mises-Fisher distributions.
 
 
 <section>
 <img src="http://visionlab.ut.ac.ir/mixest/docs/examples/img/example4.gif">
 </section>
+
 An example visualizing the estimation process of a modified version of the Competitive EM (CEM) algorithm on sample 2-D data.
 
 -------------------------------------------------------------------------------

mixest/auxillary/checking/check_grad.m

@@ -21,7 +21,7 @@
 function check_grad(D, theta, data)
 
 if nargin < 3
-    nrd = 200000;
+    nrd = 10000;
     data = D.sample(theta, nrd);
 else
     data = mxe_readdata(data);
@@ -61,7 +61,8 @@ function check_grad(D, theta, data)
 
 px = vec2obj(theta, x);
 if nargout>1
-    [y, dd] = D.ll(px, data);
+    [y, store] = D.ll(px, data);
+    dd = D.llgrad(px, data, store);
     dy = obj2vec(dd);
 else
     y = D.ll(px, data);

mixest/auxillary/manopt_manifolds/simplex/productsimplexfactory.m

@@ -0,0 +1,115 @@
+%% |simplexfactory|
+% Returns a manifold struct to optimize over simplex with given dimension
+%
+% *Syntax*
+%
+%   M = simplexfactory(n)
+%
+% *Description*
+%
+% |M = simplexfactory(n)| returns |M|, a structure describing the simplex
+% manifold of dimension |n|.
+%
+
+% Copyright 2015 Reshad Hosseini and Mohamadreza Mash'al
+% This file is part of MixEst: visionlab.ut.ac.ir/mixest
+%
+% Original author: Reshad Hosseini, Oct. 29, 2014.
+%
+% Change log: 
+%
+
+function M = productsimplexfactory(n, m)
+
+    if nargin < 2
+        m = 1;
+    end
+
+    if ~exist('n', 'var') || isempty(n)
+        n = 1;
+    end
+
+    M.name = @() sprintf('Simplex space Embedded in R^%i', n);
+
+    M.dim = @() (n-1)* m ;
+
+    M.inner = @(x, d1, d2) d1(:).'*d2(:);
+
+    M.norm = @(x, d) norm(d, 'fro');
+
+    M.dist = @(x, y) norm(x-y, 'fro');
+
+    M.typicaldist = @() sqrt(n-1);
+
+    M.proj = @(x, d) d;
+
+    M.egrad2rgrad = @egrad2rgrad;
+    function gn = egrad2rgrad(x, g)
+        gn = zeros(n-1, m);
+        for k = 1:m
+            Cov = diag(x(:,k)) - x(:,k) * x(:,k)';
+            gn(:,k) = Cov(1:end-1,:)*g(:,k);
+        end
+    end
+
+    %M.ehess2rhess = @(x, eg, eh, d) eh;
+
+    M.tangent = M.proj;
+
+    M.exp = @expmap;
+    function y = expmap(x, d, t)
+        if size(x, 1) == 1
+            y = ones(n, m);
+            return;
+        end
+        % apply first the following change of variable
+        xn = bsxfun(@minus, log(x(1:n-1, :)), log(x(n, :)));
+        % moving in the variable change domain
+        if nargin == 3
+            yn = xn + t*d;
+        else
+            yn = xn + d;
+        end
+        % going back to the original domain
+        y = exp([yn;zeros(1,m)]);
+        y = bsxfun(@rdivide, y, sum(y));
+    end
+
+    M.retr = M.exp;
+
+    M.hash = @(x) ['z' hashmd5(x(:))];
+
+    M.rand = @random;
+    function x = random()
+        x = rand(n,m);
+        x = bsxfun(@rdivide, x, sum(x));
+    end
+
+    M.randvec = @randvec;
+    function u = randvec(x) %#ok<INUSD>
+        u = randn(n-1,m);
+        u = bsxfun( @rdivide, u, sqrt(sum(u.^2)) );
+    end
+
+    M.lincomb = @lincomb;
+    function v = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
+        if nargin == 3
+            v = a1*d1;
+        elseif nargin == 5
+            v = a1*d1 + a2*d2;
+        else
+            error('Bad usage of simplex.lincomb');
+        end
+    end
+
+    M.zerovec = @(x) zeros(n-1, 1);
+
+    M.transp = @(x1, x2, d) d;
+
+    M.pairmean = @(x1, x2) .5*(x1+x2);
+
+    M.vec = @(x, u_mat) u_mat(:);
+    M.mat = @(x, u_vec) reshape(u_vec, [m, n]);
+    M.vecmatareisometries = @() true;
+
+end

mixest/auxillary/manopt_manifolds/triufactory.m

@@ -0,0 +1,85 @@
+function M = triufactory(m, n)
+% Returns a manifold struct over m-by-n lower-triangular matrices.
+%
+% function M = triufactory(m, n)
+%
+% Returns M, a structure describing the Euclidean space of lower-triangular
+% m-by-n matrices equipped with the standard Frobenius distance and 
+% associated trace inner product as a manifold for Manopt.
+
+% Copyright 2015 Reshad Hosseini and Mohamadreza Mash'al
+% This file is part of MixEst: visionlab.ut.ac.ir/mixest
+%
+% Original author: Reshad Hosseini, Aug, 08, 2015
+
+
+    if ~exist('n', 'var') || isempty(n)
+        n = 1;
+    end
+
+    M.name = @() sprintf('Euclidean space R^(%dx%d)', m, n);
+
+    M.dim = @() m*n;
+
+    M.inner = @(x, d1, d2) d1(:).'*d2(:);
+
+    M.norm = @(x, d) norm(d, 'fro');
+
+    M.dist = @(x, y) norm(x-y, 'fro');
+
+    M.typicaldist = @() sqrt(m*n);
+
+    M.proj = @(x, d) triu(d);
+
+    M.egrad2rgrad = @(x, g) triu(g);
+
+    M.ehess2rhess = @(x, eg, eh, d) eh;
+
+    M.tangent = M.proj;
+
+    M.exp = @exp;
+    function y = exp(x, d, t)
+        if nargin == 3
+            y = x + t*d;
+        else
+            y = x + d;
+        end
+        y = triu(y);
+    end
+
+    M.retr = M.exp;
+
+	M.log = @(x, y) triu(y-x);
+
+    M.hash = @(x) ['z' hashmd5(x(:))];
+
+    M.rand = @() triu(randn(m, n));
+
+    M.randvec = @randvec;
+    function u = randvec(x) %#ok<INUSD>
+        u = triu(randn(m, n));
+        u = u / norm(u, 'fro');
+    end
+
+    M.lincomb = @lincomb;
+    function v = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
+        if nargin == 3
+            v = a1*d1;
+        elseif nargin == 5
+            v = a1*d1 + a2*d2;
+        else
+            error('Bad usage of euclidean.lincomb');
+        end
+    end
+
+    M.zerovec = @(x) zeros(m, n);
+
+    M.transp = @(x1, x2, d) d;
+
+    M.pairmean = @(x1, x2) triu(.5*(x1+x2));
+
+    M.vec = @(x, u_mat) u_mat(:);
+    M.mat = @(x, u_vec) reshape(u_vec, [m, n]);
+    M.vecmatareisometries = @() true;
+
+end

mixest/distributions/mixture/mixturefactory.m

@@ -1185,7 +1185,7 @@ function fill_cache(all_data, datapatchsize)
         if ~isempty(weight)
             component_weights = bsxfun(@times, component_weights, weight);
         end
-
+        low_memory = true;
         if low_memory
             store = rmfield(store,'hX');
             %store = rmfield(store,'llik');