#NEWLINE# def min_faktor( Mu , Sigma , k , gamma = 1 ) : #NEWLINE# #INDENT# D = Mu . shape [ 0 ] #NEWLINE# logS = np . zeros ( ( ( D - 1 ) , ) ) #NEWLINE# MP = np . zeros ( ( ( D - 1 ) , ) ) #NEWLINE# P = np . zeros ( ( ( D - 1 ) , ) ) #NEWLINE# M = np . copy ( Mu ) #NEWLINE# V = np . copy ( Sigma ) #NEWLINE# b = False #NEWLINE# for count in range ( 50 ) : #NEWLINE# #INDENT# diff = 0 #NEWLINE# for i in range ( ( D - 1 ) ) : #NEWLINE# #INDENT# l = ( i if ( i < k ) else ( i + 1 ) ) #NEWLINE# try : #NEWLINE# #INDENT# ( M , V , P [ i ] , MP [ i ] , logS [ i ] , d ) = lt_factor ( k , l , M , V , MP [ i ] , P [ i ] , gamma ) #NEWLINE# #UNINDENT# except Exception as e : #NEWLINE# #INDENT# raise #NEWLINE# #UNINDENT# if np . isnan ( d ) : #NEWLINE# #INDENT# break #NEWLINE# #UNINDENT# diff += np . abs ( d ) #NEWLINE# #UNINDENT# if np . isnan ( d ) : #NEWLINE# #INDENT# break #NEWLINE# #UNINDENT# if ( np . abs ( diff ) < 0.001 ) : #NEWLINE# #INDENT# b = True #NEWLINE# break #NEWLINE# #UNINDENT# #UNINDENT# if np . isnan ( d ) : #NEWLINE# #INDENT# logZ = ( - np . Infinity ) #NEWLINE# ( yield logZ ) #NEWLINE# dlogZdMu = np . zeros ( ( D , 1 ) ) #NEWLINE# ( yield dlogZdMu ) #NEWLINE# dlogZdMudMu = np . zeros ( ( D , D ) ) #NEWLINE# ( yield dlogZdMudMu ) #NEWLINE# dlogZdSigma = np . zeros ( ( ( 0.5 * ( D * ( D + 1 ) ) ) , 1 ) ) #NEWLINE# ( yield dlogZdSigma ) #NEWLINE# mvmin = [ Mu [ k ] , Sigma [ ( k , k ) ] ] #NEWLINE# ( yield mvmin ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# C = ( np . eye ( D ) / sq2 ) #NEWLINE# C [ k , : ] = ( ( - 1 ) / sq2 ) #NEWLINE# C = np . delete ( C , k , 1 ) #NEWLINE# R = ( np . sqrt ( P . T ) * C ) #NEWLINE# r = np . sum ( ( MP . T * C ) , 1 ) #NEWLINE# mp_not_zero = np . where ( ( MP != 0 ) ) #NEWLINE# mpm = ( ( MP [ mp_not_zero ] * MP [ mp_not_zero ] ) / P [ mp_not_zero ] ) #NEWLINE# mpm = sum ( mpm ) #NEWLINE# s = sum ( logS ) #NEWLINE# IRSR = ( np . eye ( ( D - 1 ) ) + np . dot ( np . dot ( R . T , Sigma ) , R ) ) #NEWLINE# rSr = np . dot ( np . dot ( r . T , Sigma ) , r ) #NEWLINE# A = np . dot ( R , np . linalg . solve ( IRSR , R . T ) ) #NEWLINE# A = ( 0.5 * ( A . T + A ) ) #NEWLINE# b = ( Mu + np . dot ( Sigma , r ) ) #NEWLINE# Ab = np . dot ( A , b ) #NEWLINE# try : #NEWLINE# #INDENT# cIRSR = np . linalg . cholesky ( IRSR ) #NEWLINE# #UNINDENT# except np . linalg . LinAlgError : #NEWLINE# #INDENT# try : #NEWLINE# #INDENT# cIRSR = np . linalg . cholesky ( ( IRSR + ( 1e-10 * np . eye ( IRSR . shape [ 0 ] ) ) ) ) #NEWLINE# #UNINDENT# except np . linalg . LinAlgError : #NEWLINE# #INDENT# cIRSR = np . linalg . cholesky ( ( IRSR + ( 1e-06 * np . eye ( IRSR . shape [ 0 ] ) ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# dts = ( 2 * np . sum ( np . log ( np . diagonal ( cIRSR ) ) ) ) #NEWLINE# logZ = ( ( ( ( 0.5 * ( ( rSr - np . dot ( b . T , Ab ) ) - dts ) ) + np . dot ( Mu . T , r ) ) + s ) - ( 0.5 * mpm ) ) #NEWLINE# ( yield logZ ) #NEWLINE# btA = np . dot ( b . T , A ) #NEWLINE# dlogZdMu = ( r - Ab ) #NEWLINE# ( yield dlogZdMu ) #NEWLINE# dlogZdMudMu = ( - A ) #NEWLINE# ( yield dlogZdMudMu ) #NEWLINE# dlogZdSigma = ( ( ( ( - A ) - ( 2 * np . outer ( r , Ab . T ) ) ) + np . outer ( r , r . T ) ) + np . outer ( btA . T , Ab . T ) ) #NEWLINE# dlogZdSigma2 = np . zeros_like ( dlogZdSigma ) #NEWLINE# np . fill_diagonal ( dlogZdSigma2 , np . diagonal ( dlogZdSigma ) ) #NEWLINE# dlogZdSigma = ( 0.5 * ( ( dlogZdSigma + dlogZdSigma . T ) - dlogZdSigma2 ) ) #NEWLINE# dlogZdSigma = np . rot90 ( dlogZdSigma , k = 2 ) [ np . triu_indices ( D ) ] [ : : ( - 1 ) ] #NEWLINE# ( yield dlogZdSigma )