diff --git a/GPy/models/state_space_main.py b/GPy/models/state_space_main.py index dd639ad0..6ed2fbeb 100644 --- a/GPy/models/state_space_main.py +++ b/GPy/models/state_space_main.py @@ -3415,10 +3415,6 @@ class ContDescrStateSpace(DescreteStateSpace): #Q_noise = Q_noise_1 - # Return - #import pdb; pdb.set_trace() - #if dt != 0: - # Q_noise = Q_noise + np.eye(Q_noise.shape[0])*1e-8 Q_noise = 0.5*(Q_noise + Q_noise.T) # Symmetrize return A, Q_noise,None, dA, dQ @@ -3624,123 +3620,3 @@ def balance_ss_model(F,L,Qc,H,Pinf,P0,dF=None,dQc=None,dPinf=None,dP0=None): # (F,L,Qc,H,Pinf,P0,dF,dQc,dPinf,dP0) return bF, bL, bQc, bH, bPinf, bP0, bdF, bdQc, bdPinf, bdP0 - -#def psd_matrix_inverse(k,Q, U=None,S=None, p_largest_cond_num=None, regularization_type=2): -# """ -# Function inverts positive definite matrix and regularizes the inverse. -# Regularization is useful when original matrix is badly conditioned. -# Function is currently used only in SparseGP code. -# -# Inputs: -# ------------------------------ -# k: int -# Iteration umber. Used for information only. Value -1 corresponds to P_inf_inv. -# -# Q: matrix -# To be inverted -# -# U,S: matrix. vector -# SVD components of Q -# -# p_largest_cond_num: float -# Largest condition value for the inverted matrix. If cond. number is smaller than that -# no regularization happen. -# -# regularization_type: 1 or 2 -# Regularization type. -# """ -# #import pdb; pdb.set_trace() -## if (k == 0) or (k == -1): # -1 - P_inf_inv computation -## import pdb; pdb.set_trace() -# -# if p_largest_cond_num is None: -# raise ValueError("psd_matrix_inverse: None p_largest_cond_num") -# -# if U is None or S is None: -# (U, S, Vh) = sp.linalg.svd( Q, full_matrices=False, compute_uv=True, overwrite_a=False, check_finite=False) -# if S[0] < (1e-4): -# #import pdb; pdb.set_trace() -# warnings.warn("""state_space_main psd_matrix_inverse: largest singular value is too small {0:e}. -# condition number is {1:e} Maybe somethigng is wrong -# """.format(S[0], S[0]/S[-1])) -# S = S + (1e-4 - S[0]) # make the S[0] at least 1e-4 -# -# current_conditional_number = S[0]/S[-1] -# if (current_conditional_number > p_largest_cond_num): -# if (regularization_type == 1): -# regularizer = S[0] / p_largest_cond_num -# # the second computation of SVD is done to compute more precisely singular -# # vectors of small singular values, since small singular values become large. -# # It is not very clear how this step is useful but test is here. -# (U, S, Vh) = sp.linalg.svd( Q + regularizer*np.eye(Q.shape[0]), -# full_matrices=False, compute_uv=True, overwrite_a=False, check_finite=False) -# -# Q_inverse_r = np.dot( U * 1.0/S , U.T ) # Assume Q_inv is positive definite -# -# # In this case, RBF kernel we get complx eigenvalues. Probably -# # for small eigenvalue corresponding eigenvectors are not very orthogonal. -# ##########Q_inverse = np.dot( Vh.T * ( 1.0/(S + regularizer)) , U.T ) -# elif (regularization_type == 2): -# -# new_border_value = np.sqrt(current_conditional_number)/2 -# if p_largest_cond_num >= new_border_value: # this type of regularization works -# regularizer = ( S[0] / p_largest_cond_num / 2.0 )**2 -# -# Q_inverse_r = np.dot( U * ( S/(S**2 + regularizer)) , U.T ) # Assume Q_inv is positive definite -# else: -# -# better_curr_cond_num = new_border_value -# warnings.warn("""state_space_main psd_matrix_inverse: reg_type = 2 can't be done completely. -# Current conditionakl number {0:e} is reduced to {1:e} by reg_type = 1""".format(current_conditional_number, better_curr_cond_num)) -# -# regularizer = S[0] / better_curr_cond_num -# # the second computation of SVD is done to compute more precisely singular -# # vectors of small singular values, since small singular values become large. -# # It is not very clear how this step is useful but test is here. -# (U, S, Vh) = sp.linalg.svd( Q + regularizer*np.eye(Q.shape[0]), -# full_matrices=False, compute_uv=True, overwrite_a=False, check_finite=False) -# -# regularizer = ( S[0] / p_largest_cond_num / 2.0 )**2 -# -# Q_inverse_r = np.dot( U * ( S/(S**2 + regularizer)) , U.T ) # Assume Q_inv is positive definite -# -# assert regularizer*10 < S[0], "regularizer is not << S[0]" -# assert regularizer > 10*S[-1], "regularizer is not >> S[-1]" -# -## Old version -> -## lamda_star = np.sqrt(current_conditional_number) -## if 2*p_largest_cond_num >= lamda_star: -## lamda = current_conditional_number / 2 / p_largest_cond_num -## -## regularizer = (S[-1] * lamda)**2 -## -## Q_inverse_r = np.dot( U * ( S/(S**2 + regularizer)) , U.T ) # Assume Q_inv is positive definite -## else: -## better_curr_cond_num = (2*p_largest_cond_num)**2 / 2 # division by 2 just in case here -## warnings.warn("""state_space_main psd_matrix_inverse: reg_type = 2 can't be done completely. -## Current conditionakl number {0:e} is reduced to {1:e} by reg_type = 1""".format(current_conditional_number, better_curr_cond_num)) -## -## regularizer = S[0] / better_curr_cond_num -## # the second computation of SVD is done to compute more precisely singular -## # vectors of small singular values, since small singular values become large. -## # It is not very clear how this step is useful but test is here. -## (U, S, Vh) = sp.linalg.svd( Q + regularizer*np.eye(Q.shape[0]), -## full_matrices=False, compute_uv=True, overwrite_a=False, check_finite=False) -## -## lamda = better_curr_cond_num / 2 / p_largest_cond_num -## -## regularizer = (S[-1] * lamda)**2 -## -## Q_inverse_r = np.dot( U * ( S/(S**2 + regularizer)) , U.T ) # Assume Q_inv is positive definite -## -## assert lamda > 10, "Some assumptions are incorrect if this is not satisfied." -## Old version <- -# else: -# raise ValueError("AQcompute_batch_Python:Q_inverse: Invalid regularization type") -# -# else: -# Q_inverse_r = np.dot( U * 1.0/S , U.T ) # Assume Q_inv is positive definite -# # When checking conditional number 2 times difference is ok. -# Q_inverse_r = 0.5*(Q_inverse_r + Q_inverse_r.T) -# -# return Q_inverse_r \ No newline at end of file diff --git a/GPy/models/state_space_model.py b/GPy/models/state_space_model.py index c74a25f2..d16b5adc 100644 --- a/GPy/models/state_space_model.py +++ b/GPy/models/state_space_model.py @@ -100,8 +100,8 @@ class StateSpace(Model): # Get the model matrices from the kernel (F,L,Qc,H,P_inf, P0, dFt,dQct,dP_inft, dP0t) = self.kern.sde() - #Qc = Qc + np.eye(Qc.shape[0]) * 1e-8 - #import pdb; pdb.set_trace() + + # necessary parameters measurement_dim = self.output_dim grad_params_no = dFt.shape[2]+1 # we also add measurement noise as a parameter