<jupyter_start><jupyter_text>Latent Dirichlet Allocation- Optimization Part To optimize our model, we consider two different aspects- Algorithm Optimization and Code Optimization  1. Algorithm Optimization  Optimize the Normalization Process After updating $\phi$, we need to normalize $\phi$ to make it a probabiliy.So we have,$$\phi_i = \frac{\phi_i}{\phi1+\phi_2+...+\phi_K}$$Actually we can use $\log$ function to optimize this process,$$log(\phi_i) = \log(\phi_i)-\log(\phi_1+\phi_2+...+\phi_K$$So in normalization, this method can add in pairs during each iteration after geting $log(\phi_1), log(\phi_2),...,log(\phi_K)$ and finally speed up the summation of $log(\phi_1+\phi_2+...+\phi_K)$And because of the target value has become $log(\phi_1)$, the iterative formula of $\phi$ should be changed into:$$log(\phi_{dni}^{t+1}) = log(\beta_{i,w_n})+\psi(\gamma_{di}^t)$$<jupyter_code># for i in range(K):   
#     log_phi[m][n,i] = np.log(beta[i,np.int(words[m][n])]) + digamma(gamma[m,i])                
#     logsum = log_sum(logsum, log_phi[m][n,i])
# log_phi_mn = log_phi[m][n,:] - logsum              # use new metohd to implement nomalization
# log_phi[m][n,:] = log_phi_mn
                
# phi[m][n,:] = np.exp(log_phi_mn)<jupyter_output><empty_output><jupyter_text>Optimize $\gamma$ Updating Process When we are updating the $\gamma$, we use$$\gamma_i^0 = \alpha_i+\frac{N}{K}\\\gamma_i^{t+1} = \alpha_i+\sum_{n=1}^N \phi_{dn}^{t+1}$$But we can find that the $\gamma^t$ has the $\alpha$ information $\gamma^{t+1}$ needs, so we can optimize this process$$\begin{aligned}\gamma_i^{t+1} &= \alpha_i+\sum_{n=1}^N \phi_{dn}^{t+1}\\&=\alpha_i+\sum_{n=1}^N (\phi_n^{t}+\bigtriangleup \phi_{dn}^{t+1})\\&=(\alpha_i+\sum_{n=1}^N \phi_n^{t})+\sum_{n=1}^N\bigtriangleup \phi_{dn}^{t+1})\\&=\gamma_i^t+\sum_{n=1}^N\left(\phi_{dn}^{t+1}-\phi_{dn}^{t}\right)\end{aligned}$$<jupyter_code># d_phi = phi[m] - phi_old[m]
# gamma[m,:]  = gamma[m,:] + np.sum(d_phi, axis = 0)<jupyter_output><empty_output><jupyter_text>Vectorization We also use vector/matrix manipulation to replace the for loop<jupyter_code># Example
#
# for-loop
# for i in range(K):
#             g1 = M*(digamma(np.sum(alpha))-digamma(alpha[i]))
#             g2 = 0
#             for d in range(M):
#                 g2 += digamma(gamma[d,i])-digamma(np.sum(gamma[d,:]))
#             g[i] = g1 + g2
# Vector
# g = M*(digamma(np.sum(alpha))-digamma(alpha)) 
#         + np.sum(digamma(gamma) -np.tile(digamma(np.sum(gamma,axis=1)),(K,1)).T,axis=0)<jupyter_output><empty_output><jupyter_text>The Complete Code for Algorithm Optimization Process<jupyter_code># a new function to calculate log of sum
def log_sum(log_a, log_b):
    """
    input: log(a), log(b)
    output: log(a+b)
    """
    return log_a + np.log(1+np.exp(log_b - log_a))



def optimize_vp_opt(phi, gamma, alpha, beta, words, M, N, K, max_iter=500):
    '''
    optimize the variational parameter
    
    Parameters
    ----------
    phi:   ndarray
           An array of topic-word matrix
    gamma: ndarray
           A matrix of doc-topic
    alpha: ndarray
           the parameter of doc-topic dirichlet distribution
    beta:  ndarray
           the parameter of topic-word dirichlet distribution
    words: list 
           the list of lists of words in all 
    M : int, the number of documents
    N : ndarraay, the number of words in each document
    K : int, the number of topics in the corpus
    Returns
    -------
    out : list of ndarray
          the optimized and normalized(sum to 1) phi 
    '''
    
    for t in range(max_iter):
        phi_old = phi
        
        # we use log(phi) here and following processes
        log_phi = np.array(list(map(np.log, phi)))
        gamma_old = gamma
       
        for m in range(M):
            for n in range(N[m]):
                
                logsum = 0
                for i in range(K):
                    
                    # use new method in log form to update phi
                    log_phi[m][n,i] = np.log(beta[i,np.int(words[m][n])]) + digamma(gamma[m,i])
                    
                    logsum = log_sum(logsum, log_phi[m][n,i])
                # use new metohd to implement nomalization
                log_phi_mn = log_phi[m][n,:] - logsum
                log_phi[m][n,:] = log_phi_mn
                
                phi[m][n,:] = np.exp(log_phi_mn)
        
            # instead of alpha, use old phi and new phi to iterative
            d_phi = phi[m] - phi_old[m]
            gamma[m,:]  = gamma[m,:] + np.sum(d_phi, axis = 0)
            
        phi_new = phi
        gamma_new = gamma
        
        if is_convergence1(phi_old, phi_new) == True and is_convergence2(gamma_old, gamma_new) == True:
            break
   
    return phi, gamma
# estimate alpha
def alpha_estimate_opt(gamma, alpha_initial, K, M, max_iter = 100):
    """
    This is an estimation function, especially used in the process of LDA algorithm.
    digamma function and polygamma function are used in the following process.
    
    input:
    alpha_initial: the initial setting of alpha, it is an 1*K vector
    K: the number of topics
    M: the number of documents
    gamma: the result from another update function (see gamma_update())
    """
    from scipy.special import digamma, polygamma
    
    alpha = alpha_initial
    for t in range(max_iter):
        alpha_old = alpha
        
        # we use vector instead of calculating in loop
        g = M*(digamma(np.sum(alpha))-digamma(alpha)) 
        + np.sum(digamma(gamma) -np.tile(digamma(np.sum(gamma,axis=1)),(K,1)).T,axis=0)
        h = -M*polygamma(1,alpha)
        
        z = M*polygamma(1, np.sum(alpha))
        c = (np.sum(g/h))/(z**(-1) + np.sum(h**(-1)))
                           
        # update alpha                   
        alpha -= (g-c)/h
        
        if is_convergence2(alpha_old, alpha):
            break
            
    return alpha

# estimate beta
def beta_estimate_opt(K, V_words, phi, D):
    
    """
    This is an estimation function, especially used in the process of LDA algorithm
    
    input:
    K: the number of topics
    V_words: a vector of all unique words in the vocabulary
    D: D = (w_1,w_2,...w_M), contains all words in all documents
    phi: the result from another update function (see phi_update())
    
    output:
    beta: the estimate parameter for LDA, it is a K*V matrix
    """
    V = len(V_words)
    beta = np.ones((K,V))
    # first obtain the propotion values
    for j in range(V):
        word = V_words[j]
        # give a TRUE or FALSE "matrix", remember w_mnj should have the same shape with phi
        w_mnj = [np.repeat(w==word, K).reshape((len(w),K)) for w in D]
        # compute the inner sum over number of words
        sum1 = list(map(lambda x: np.sum(x,axis=0),phi*w_mnj))
        # compute the outer sum over documents
        beta[:,j] = np.sum(np.array(sum1), axis = 0)
    
    # then normalize each row s.t. the row sum is one, in vector method
    beta= beta/ np.sum(beta, axis = 1).reshape((-1,1))
        
    return beta<jupyter_output><empty_output><jupyter_text>Test<jupyter_code>import lda_package     #our package
import numpy as np
#Preprocess data
doc_a = "Brocolli is good to eat. My brother likes to eat good brocolli, but not my mother."
doc_b = "My mother spends a lot of time driving my brother around to baseball practice."
doc_c = "Some health experts suggest that driving may cause increased tension and blood pressure."
doc_d = "I often feel pressure to perform well at school, but my mother never seems to drive my brother to do better."
doc_e = "Health professionals say that brocolli is good for your health." 
doc_set = [doc_a, doc_b, doc_c, doc_d, doc_e]
texts, dictionary, corpus = lda_package.data_clean(doc_set)
text_ = lda_package.data_process(texts, dictionary)
#initialize parameter
M = len(texts)
k = 2
N = np.array(list(map(len, text_)))
V = len(dictionary)
V_words = range(V)
alpha = np.random.dirichlet(10*np.ones(k),1)[0]
beta = np.random.dirichlet(np.ones(V),k)

phi = np.array([1/k*np.ones([N[m],k]) for m in range(M)])
gamma = np.tile(alpha,(M,1)) + np.tile(N/k,(k,1)).T
%%time
ans1 = lda_package.variation_EM(M, k, text_, N, V_words, alpha, beta, gamma, phi, iteration = 1000)
%%time
ans2 = lda_package.variation_EM_new(M, k, text_, N, V_words, alpha, beta, gamma, phi, iteration = 1000)<jupyter_output>CPU times: user 12 ms, sys: 0 ns, total: 12 ms
Wall time: 10.8 ms<jupyter_text>Parallelism Optimization in numba<jupyter_code>from numba import jit, prange
from scipy.special import digamma, polygamma
@jit(nopython=True, nogil=True)
def log_sum_jit(log_a, log_b):
    """
    input: log(a), log(b)
    output: log(a+b)
    """
    return log_a + np.log(1+np.exp(log_b - log_a))

@jit(nopython=True, parallel=True)
def optimize_vp_opt_jit(phi, gamma, alpha, beta, words, M, N, K, max_iter=500):
    '''
    optimize the variational parameter
    
    Parameters
    ----------
    phi:   ndarray
           An array of topic-word matrix
    gamma: ndarray
           A matrix of doc-topic
    alpha: ndarray
           the parameter of doc-topic dirichlet distribution
    beta:  ndarray
           the parameter of topic-word dirichlet distribution
    words: list 
           the list of lists of words in all 
    M : int, the number of documents
    N : ndarraay, the number of words in each document
    K : int, the number of topics in the corpus
    Returns
    -------
    out : list of ndarray
          the optimized and normalized(sum to 1) phi 
    '''
    
    for t in range(max_iter):
        phi_old = phi
        
        # we use log(phi) here and following processes
        log_phi = np.array(list(map(np.log, phi)))
        gamma_old = gamma
       
        for m in prange(M):
            for n in prange(N[m]):
                
                logsum = 0
                for i in prange(K):
                    
                    # use new method in log form to update phi
                    log_phi[m][n,i] = np.log(beta[i,np.int(words[m][n])]) + digamma(gamma[m,i])
                    
                    logsum = log_sum(logsum, log_phi[m][n,i])
                # use new metohd to implement nomalization
                log_phi_mn = log_phi[m][n,:] - logsum
                log_phi[m][n,:] = log_phi_mn
                
                phi[m][n,:] = np.exp(log_phi_mn)
        
            # instead of alpha, use old phi and new phi to iterative
            d_phi = phi[m] - phi_old[m]
            gamma[m,:]  = gamma[m,:] + np.sum(d_phi, axis = 0)
            
        phi_new = phi
        gamma_new = gamma
        
        if is_convergence1_jit(phi_old, phi_new) == True and is_convergence2_jit(gamma_old, gamma_new) == True:
            break
   
    return phi, gamma


def is_convergence1_jit(old, new, tol = 10**(-2)):
    """
    output:
    TRUR or FALSE
    """
    loss = np.sqrt(list(map(np.sum,np.square(old - new))))
    return np.max(loss) <= tol

def is_convergence2_jit(old, new, tol = 10**(-2)):
    """
    output:
    TRUR or FALSE
    """
    loss = np.sqrt(np.sum(np.square(old - new)))
    return np.max(loss) <= tol

# estimate alpha
def alpha_estimate_opt_jit(gamma, alpha_initial, K, M, max_iter = 100):
    """
    This is an estimation function, especially used in the process of LDA algorithm.
    digamma function and polygamma function are used in the following process.
    
    input:
    alpha_initial: the initial setting of alpha, it is an 1*K vector
    K: the number of topics
    M: the number of documents
    gamma: the result from another update function (see gamma_update())
    """
    from scipy.special import digamma, polygamma
    
    alpha = alpha_initial
    for t in range(max_iter):
        alpha_old = alpha
        
        # we use vector instead of calculating in loop
        g = M*(digamma(np.sum(alpha))-digamma(alpha)) 
        + np.sum(digamma(gamma) -np.tile(digamma(np.sum(gamma,axis=1)),(K,1)).T,axis=0)
        h = -M*polygamma(1,alpha)
        
        z = M*polygamma(1, np.sum(alpha))
        c = (np.sum(g/h))/(z**(-1) + np.sum(h**(-1)))
                           
        # update alpha                   
        alpha -= (g-c)/h
        
        if is_convergence2_jit(alpha_old, alpha):
            break
            
    return alpha


# estimate beta
def beta_estimate_opt_jit(K, V_words, phi, D):
    
    """
    This is an estimation function, especially used in the process of LDA algorithm
    
    input:
    K: the number of topics
    V_words: a vector of all unique words in the vocabulary
    D: D = (w_1,w_2,...w_M), contains all words in all documents
    phi: the result from another update function (see phi_update())
    
    output:
    beta: the estimate parameter for LDA, it is a K*V matrix
    """
    V = len(V_words)
    beta = np.ones((K,V))
    # first obtain the propotion values
    for j in prange(V):
        word = V_words[j]
        # give a TRUE or FALSE "matrix", remember w_mnj should have the same shape with phi
        w_mnj = [np.repeat(w==word, K).reshape((len(w),K)) for w in D]
        # compute the inner sum over number of words
        sum1 = list(map(lambda x: np.sum(x,axis=0),phi*w_mnj))
        # compute the outer sum over documents
        beta[:,j] = np.sum(np.array(sum1), axis = 0)
    
    # then normalize each row s.t. the row sum is one, in vector method
    beta= beta/ np.sum(beta, axis = 1).reshape((-1,1))
        
    return beta


def variation_EM_new_jit(M, K, D, N, V_words, alpha_initial, beta_initial, gamma_initial, phi_initial, iteration = 1000):
    
    phi_gamma = optimize_vp_opt_jit(phi_initial, gamma_initial, alpha_initial, beta_initial, D, M, N, K)
    phi = phi_gamma[0]
    gamma = phi_gamma[1]
    
     
    (alpha, beta) = (alpha_initial, beta_initial)
    
    for t in range(iteration):
        
        (phi_old, gamma_old) = (phi, gamma)
        
        alpha = alpha_estimate_opt_jit(gamma, alpha, K, M)
        beta = beta_estimate_opt_jit(K, V_words, phi, D)
        
        phi_gamma1 = optimize_vp_opt_jit(phi, gamma, alpha, beta, D, M, N, K)
        phi = phi_gamma1[0]
        gamma = phi_gamma1[1]
        
        if is_convergence2_jit(gamma_old, gamma) and is_convergence1_jit(phi_old, phi):
            break
    
    return alpha, beta, gamma, phi
%%time
ans3 = lda_package.variation_EM_new(M, k, text_, N, V_words, alpha, beta, gamma, phi, iteration = 1000)
%%time
ans4 = variation_EM_new_jit(M, k, text_, N, V_words, alpha, beta, gamma, phi, iteration = 1000)<jupyter_output>CPU times: user 8 ms, sys: 4 ms, total: 12 ms
Wall time: 8.84 ms
{'fraction_non_alphanumeric': 0.12337263990184719, 'fraction_numerical': 0.010905868720605276, 'mean_word_length': 2.449800141076887, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '278321', 'license': 'MIT', 'n_tokens_mistral': 5168, 'n_tokens_neox': 4921, 'n_words': 1433, 'path': 'LDA-Optimization.ipynb', 'repo_name': 'yangbaovera/663-Final-Project', 'size': '14671'}