Computes the cholesky decomposition of a symetric semidefinite matrix or of a definite positive matrix.
0001 function T = reduced_rank_cholesky(X) 0002 % Computes the cholesky decomposition of a symetric semidefinite matrix or of a definite positive matrix. 0003 0004 %@info: 0005 %! @deftypefn {Function File} { @var{T} =} reduced_rank_cholesky (@var{X}) 0006 %! @anchor{reduced_rank_cholesky} 0007 %! @sp 1 0008 %! Computes the cholesky decomposition of a symetric semidefinite matrix or of a definite positive matrix. 0009 %! @sp 2 0010 %! @strong{Inputs} 0011 %! @sp 1 0012 %! @table @ @var 0013 %! @item X 0014 %! n*n matrix of doubles to be factorized (X is supposed to be semidefinite positive). 0015 %! @end table 0016 %! @sp 2 0017 %! @strong{Outputs} 0018 %! @sp 1 0019 %! @table @ @var 0020 %! @item T 0021 %! q*n matrix of doubles such that T'*T = X, where q is the number of positive eigenvalues in X. 0022 %! @end table 0023 %! @sp 2 0024 %! @strong{Remarks} 0025 %! @sp 1 0026 %! [1] If X is not positive definite, then X has to be a symetric semidefinite matrix. 0027 %! @sp 1 0028 %! [2] The matrix T is upper triangular iff X is positive definite. 0029 %! @sp 2 0030 %! @strong{This function is called by:} 0031 %! @sp 1 0032 %! @ref{particle/sequential_importance_particle_filter} 0033 %! @sp 2 0034 %! @strong{This function calls:} 0035 %! @sp 2 0036 %! @end deftypefn 0037 %@eod: 0038 0039 % Copyright (C) 2009, 2010, 2011 Dynare Team 0040 % 0041 % This file is part of Dynare. 0042 % 0043 % Dynare is free software: you can redistribute it and/or modify 0044 % it under the terms of the GNU General Public License as published by 0045 % the Free Software Foundation, either version 3 of the License, or 0046 % (at your option) any later version. 0047 % 0048 % Dynare is distributed in the hope that it will be useful, 0049 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0050 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0051 % GNU General Public License for more details. 0052 % 0053 % You should have received a copy of the GNU General Public License 0054 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0055 0056 [T,X_is_not_positive_definite] = chol(X); 0057 0058 if X_is_not_positive_definite 0059 n = length(X); 0060 [U,D] = eig(X); 0061 [tmp,max_elements_indices] = max(abs(U),[],1); 0062 negloc = (U(max_elements_indices+(0:n:(n-1)*n))<0); 0063 U(:,negloc) = -U(:,negloc); 0064 D = diag(D); 0065 tol = sqrt(eps(max(D))*length(D)*10); 0066 t = (abs(D) > tol); 0067 D = D(t); 0068 if ~(sum(D<0)) 0069 T = diag(sqrt(D))*U(:,t)'; 0070 else 0071 disp('reduced_rank_cholesky:: Input matrix is not semidefinite positive!') 0072 T = NaN; 0073 end 0074 end 0075 0076 %@test:1 0077 %$ n = 10; 0078 %$ m = 100; 0079 %$ 0080 %$ X = randn(n,m); 0081 %$ X = X*X'; 0082 %$ 0083 %$ t = ones(2,1); 0084 %$ 0085 %$ try 0086 %$ T = reduced_rank_cholesky(X); 0087 %$ catch 0088 %$ t(1) = 0; 0089 %$ T = all(t); 0090 %$ return; 0091 %$ end 0092 %$ 0093 %$ 0094 %$ % Check the results. 0095 %$ t(2) = dyn_assert(T,chol(X),1e-16); 0096 %$ T = all(t); 0097 %@eof:1