def EVBMF(Y, sigma2=None, H=None):
    'Implementation of the analytical solution to Empirical Variational\n        Bayes Matrix Factorization.\n\n    This function can be used to calculate the analytical solution to\n        empirical VBMF.\n    This is based on the paper and MatLab code by Nakajima et al.:\n    "Global analytic solution of fully-observed variational Bayesian matrix\n        factorization."\n\n    Notes\n    -----\n        If sigma2 is unspecified, it is estimated by minimizing the free\n            energy.\n        If H is unspecified, it is set to the smallest of the sides of the\n            input Y.\n\n    Attributes\n    ----------\n    Y : numpy-array\n        Input matrix that is to be factorized. Y has shape (L,M), where L<=M.\n\n    sigma2 : int or None (default=None)\n        Variance of the noise on Y.\n\n    H : int or None (default = None)\n        Maximum rank of the factorized matrices.\n\n    Returns\n    -------\n    U : numpy-array\n        Left-singular vectors.\n\n    S : numpy-array\n        Diagonal matrix of singular values.\n\n    V : numpy-array\n        Right-singular vectors.\n\n    post : dictionary\n        Dictionary containing the computed posterior values.\n\n\n    References\n    ----------\n    .. [1] Nakajima, Shinichi, et al. "Global analytic solution of\n        fully-observed variational Bayesian matrix factorization." Journal of\n        Machine Learning Research 14.Jan (2013): 1-37.\n\n    .. [2] Nakajima, Shinichi, et al. "Perfect dimensionality recovery by\n        variational Bayesian PCA." Advances in Neural Information Processing\n        Systems. 2012.\n    '
    (L, M) = Y.shape
    if (H is None):
        H = L
    alpha = (L / M)
    tauubar = (2.5129 * np.sqrt(alpha))
    (U, s, V) = torch.svd(Y)
    U = U[:, :H]
    s = s[:H]
    V = V[:H].T
    residual = 0.0
    if (H < L):
        residual = torch.sum((np.sum((Y ** 2)) - np.sum((s ** 2))))
    if (sigma2 is None):
        xubar = ((1 + tauubar) * (1 + (alpha / tauubar)))
        eH_ub = (int(np.min([(np.ceil((L / (1 + alpha))) - 1), H])) - 1)
        upper_bound = ((torch.sum((s ** 2)) + residual) / (L * M))
        lower_bound = torch.max(torch.stack([((s[(eH_ub + 1)] ** 2) / (M * xubar)), (torch.mean((s[(eH_ub + 1):] ** 2)) / M)], dim=0))
        scale = 1.0
        s = (s * np.sqrt(scale))
        residual = (residual * scale)
        lower_bound = (lower_bound * scale)
        upper_bound = (upper_bound * scale)
        sigma2_opt = minimize_scalar(EVBsigma2, args=(L, M, s.cpu().numpy(), residual, xubar), bounds=[lower_bound.cpu().numpy(), upper_bound.cpu().numpy()], method='Bounded')
        sigma2 = sigma2_opt.x
    threshold = np.sqrt((((M * sigma2) * (1 + tauubar)) * (1 + (alpha / tauubar))))
    pos = torch.sum((s > threshold))
    d = ((s[:pos] / 2) * ((1 - (((L + M) * sigma2) / (s[:pos] ** 2))) + torch.sqrt((((1 - (((L + M) * sigma2) / (s[:pos] ** 2))) ** 2) - ((((4 * L) * M) * (sigma2 ** 2)) / (s[:pos] ** 4))))))
    return (U[:, :pos], torch.diag(d), V[:, :pos])
Implementation of the analytical solution to Empirical Variational
    Bayes Matrix Factorization.

This function can be used to calculate the analytical solution to
    empirical VBMF.
This is based on the paper and MatLab code by Nakajima et al.:
"Global analytic solution of fully-observed variational Bayesian matrix
    factorization."

Notes
-----
    If sigma2 is unspecified, it is estimated by minimizing the free
        energy.
    If H is unspecified, it is set to the smallest of the sides of the
        input Y.

Attributes
----------
Y : numpy-array
    Input matrix that is to be factorized. Y has shape (L,M), where L<=M.

sigma2 : int or None (default=None)
    Variance of the noise on Y.

H : int or None (default = None)
    Maximum rank of the factorized matrices.

Returns
-------
U : numpy-array
    Left-singular vectors.

S : numpy-array
    Diagonal matrix of singular values.

V : numpy-array
    Right-singular vectors.

post : dictionary
    Dictionary containing the computed posterior values.


References
----------
.. [1] Nakajima, Shinichi, et al. "Global analytic solution of
    fully-observed variational Bayesian matrix factorization." Journal of
    Machine Learning Research 14.Jan (2013): 1-37.

.. [2] Nakajima, Shinichi, et al. "Perfect dimensionality recovery by
    variational Bayesian PCA." Advances in Neural Information Processing
    Systems. 2012.
def EVBMF(Y, sigma2=None, H=None):
    'Implementation of the analytical solution to Empirical Variational\n        Bayes Matrix Factorization.\n\n    This function can be used to calculate the analytical solution to\n        empirical VBMF.\n    This is based on the paper and MatLab code by Nakajima et al.:\n    "Global analytic solution of fully-observed variational Bayesian matrix\n        factorization."\n\n    Notes\n    -----\n        If sigma2 is unspecified, it is estimated by minimizing the free\n            energy.\n        If H is unspecified, it is set to the smallest of the sides of the\n            input Y.\n\n    Attributes\n    ----------\n    Y : numpy-array\n        Input matrix that is to be factorized. Y has shape (L,M), where L<=M.\n\n    sigma2 : int or None (default=None)\n        Variance of the noise on Y.\n\n    H : int or None (default = None)\n        Maximum rank of the factorized matrices.\n\n    Returns\n    -------\n    U : numpy-array\n        Left-singular vectors.\n\n    S : numpy-array\n        Diagonal matrix of singular values.\n\n    V : numpy-array\n        Right-singular vectors.\n\n    post : dictionary\n        Dictionary containing the computed posterior values.\n\n\n    References\n    ----------\n    .. [1] Nakajima, Shinichi, et al. "Global analytic solution of\n        fully-observed variational Bayesian matrix factorization." Journal of\n        Machine Learning Research 14.Jan (2013): 1-37.\n\n    .. [2] Nakajima, Shinichi, et al. "Perfect dimensionality recovery by\n        variational Bayesian PCA." Advances in Neural Information Processing\n        Systems. 2012.\n    '
    (L, M) = Y.shape
    if (H is None):
        H = L
    alpha = (L / M)
    tauubar = (2.5129 * np.sqrt(alpha))
    (U, s, V) = torch.svd(Y)
    U = U[:, :H]
    s = s[:H]
    V = V[:H].T
    residual = 0.0
    if (H < L):
        residual = torch.sum((np.sum((Y ** 2)) - np.sum((s ** 2))))
    if (sigma2 is None):
        xubar = ((1 + tauubar) * (1 + (alpha / tauubar)))
        eH_ub = (int(np.min([(np.ceil((L / (1 + alpha))) - 1), H])) - 1)
        upper_bound = ((torch.sum((s ** 2)) + residual) / (L * M))
        lower_bound = torch.max(torch.stack([((s[(eH_ub + 1)] ** 2) / (M * xubar)), (torch.mean((s[(eH_ub + 1):] ** 2)) / M)], dim=0))
        scale = 1.0
        s = (s * np.sqrt(scale))
        residual = (residual * scale)
        lower_bound = (lower_bound * scale)
        upper_bound = (upper_bound * scale)
        sigma2_opt = minimize_scalar(EVBsigma2, args=(L, M, s.cpu().numpy(), residual, xubar), bounds=[lower_bound.cpu().numpy(), upper_bound.cpu().numpy()], method='Bounded')
        sigma2 = sigma2_opt.x
    threshold = np.sqrt((((M * sigma2) * (1 + tauubar)) * (1 + (alpha / tauubar))))
    pos = torch.sum((s > threshold))
    d = ((s[:pos] / 2) * ((1 - (((L + M) * sigma2) / (s[:pos] ** 2))) + torch.sqrt((((1 - (((L + M) * sigma2) / (s[:pos] ** 2))) ** 2) - ((((4 * L) * M) * (sigma2 ** 2)) / (s[:pos] ** 4))))))
    return (U[:, :pos], torch.diag(d), V[:, :pos])