1function T = reduced_rank_cholesky(X) 2% Computes the cholesky decomposition of a symetric semidefinite matrix or of a definite positive matrix. 3 4%@info: 5%! @deftypefn {Function File} { @var{T} =} reduced_rank_cholesky (@var{X}) 6%! @anchor{reduced_rank_cholesky} 7%! @sp 1 8%! Computes the cholesky decomposition of a symetric semidefinite matrix or of a definite positive matrix. 9%! @sp 2 10%! @strong{Inputs} 11%! @sp 1 12%! @table @ @var 13%! @item X 14%! n*n matrix of doubles to be factorized (X is supposed to be semidefinite positive). 15%! @end table 16%! @sp 2 17%! @strong{Outputs} 18%! @sp 1 19%! @table @ @var 20%! @item T 21%! q*n matrix of doubles such that T'*T = X, where q is the number of positive eigenvalues in X. 22%! @end table 23%! @sp 2 24%! @strong{Remarks} 25%! @sp 1 26%! [1] If X is not positive definite, then X has to be a symetric semidefinite matrix. 27%! @sp 1 28%! [2] The matrix T is upper triangular iff X is positive definite. 29%! @sp 2 30%! @strong{This function is called by:} 31%! @sp 1 32%! @ref{particle/sequential_importance_particle_filter} 33%! @sp 2 34%! @strong{This function calls:} 35%! @sp 2 36%! @end deftypefn 37%@eod: 38 39% Copyright (C) 2009-2017 Dynare Team 40% 41% This file is part of Dynare. 42% 43% Dynare is free software: you can redistribute it and/or modify 44% it under the terms of the GNU General Public License as published by 45% the Free Software Foundation, either version 3 of the License, or 46% (at your option) any later version. 47% 48% Dynare is distributed in the hope that it will be useful, 49% but WITHOUT ANY WARRANTY; without even the implied warranty of 50% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 51% GNU General Public License for more details. 52% 53% You should have received a copy of the GNU General Public License 54% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 55 56[T,X_is_not_positive_definite] = chol(X); 57 58if X_is_not_positive_definite 59 n = length(X); 60 [U,D] = eig(X); 61 [tmp,max_elements_indices] = max(abs(U),[],1); 62 negloc = (U(max_elements_indices+(0:n:(n-1)*n))<0); 63 U(:,negloc) = -U(:,negloc); 64 D = diag(D); 65 tol = sqrt(eps(max(D))*length(D)*10); 66 t = (abs(D) > tol); 67 D = D(t); 68 if ~(sum(D<0)) 69 T = diag(sqrt(D))*U(:,t)'; 70 else 71 disp('reduced_rank_cholesky:: Input matrix is not semidefinite positive!') 72 T = NaN; 73 end 74end 75 76%@test:1 77%$ n = 10; 78%$ m = 100; 79%$ 80%$ X = randn(n,m); 81%$ X = X*X'; 82%$ 83%$ t = ones(2,1); 84%$ 85%$ try 86%$ T = reduced_rank_cholesky(X); 87%$ catch 88%$ t(1) = 0; 89%$ T = all(t); 90%$ return 91%$ end 92%$ 93%$ 94%$ % Check the results. 95%$ t(2) = dassert(T,chol(X),1e-16); 96%$ T = all(t); 97%@eof:1