Implement `LogisticPCA` as a distinct variant of matrix decomposition useful for binary data

### Describe the workflow you want to enable

Currently, there is no included implementation of a PCA algorithm made for handling binary data in the scikit-learn library. However, the algorithm for "logistic PCA" is well founded, although different methods for its estimation exist, and it meets the inclusion criteria of being at least 3 years since first publication, having over 200 citations between known papers, and has a wide use and usefulness. The proposed implementation follows that of Lee et al. (2010), which allows for the use of an optional L1 regularization term. The intent is to create a `LogisticPCA` implementation that mimics the existing `SparsePCA` implementation to the greatest degree possible acknowledging their important differences, and which utilizes the adopted best practices to fit within the existing matrix decomposition algorithms and the current API of scikit-learn.

References

Tipping, Michael E. "Probabilistic visualisation of high-dimensional binary data." Advances in Neural Information Processing Systems (1999): 592-598.

Lee, Seokho, Jianhua Z. Huang, and Jianhua Hu. "Sparse logistic principal components analysis for binary data." The Annals of Applied Statistics, 4.3 (2010): 1579-1601.

Landgraf, Andrew J. and Yoonkyung Lee. "Dimensionality reduction for binary data through the projection of natural parameters." Journal of Multivariate Analysis, v.180 (2020): 104668.

### Describe your proposed solution

The sigmoid function, defined as `sig = lambda x: 1 / (1 + np.exp(-x))`, plays a crucial role in the logistic PCA algorithm. Its purpose is to map any real-valued number into the range (0, 1), which is particularly useful for binary data. 

Here are the roles the sigmoid function plays in the algorithm:

Probability Estimation: The sigmoid function is used to estimate probabilities. In the context of logistic PCA, it helps in modeling the probability that a given binary variable is 1. This is essential for handling binary data, where each element in the data matrix `X` is either 0 or 1.

Variational Parameters: The lambda function `lam = lambda x: (0.5 - sig(x)) / (2 * x)` uses the sigmoid function to compute a term that is part of the variational parameters. These parameters are crucial for approximating the posterior distribution of the latent variables.

Gradient and Hessian Computation: During the optimization steps, the sigmoid function is used to compute gradients and Hessians. These are necessary for updating the model parameters (`W` and `b`) to minimize the reconstruction error.

By transforming the data through the sigmoid function, the logistic PCA algorithm can effectively handle the binary nature of the input data, making it suitable for applications where the data is represented as binary variables.

Here are the key steps involved:

Initialization:

- Initialize the identity matrix `I`, the projection matrix `W`, latent variables `mu`, bias `b`, covariance matrices `C`, and variational parameters `xi`. Using the constants `n_components`, `n_samples`, and `n_features` where applicable in their definitions.

- Define the sigmoid function `sig` and the lambda function `lam`.

Iterative Optimization:

The optimization process runs for a specified number of iterations to `max_iter`.

Step 1. Update Sufficient Statistics:
- For each sample `n`:
   - Extract the sample `x_n` from the data matrix `X`.
   - Compute the approximation `lam_n` using the lambda function `lam`.
   - Update the covariance matrix `C` and the latent variables `mu` using the computed values.

Step 2. Optimize Variational Parameters:
- For each sample `n`:
   - Compute the posterior statistics `z` and `E_zz`.
   - Calculate `xixi`, which involves the projection matrix `W`, the posterior statistics, and the bias `b`.
   - Update the variational parameters `xi` based on the computed values.

Step 3. Update Model Parameters:
- Compute the expected value `E_zhzh` for each sample `n`.
- Compute the expected value `E_zh`.
- For each dimension `i`:
   - Compute the approximation `lam_i` using the lambda function `lam`.
   - Calculate the gradient `g` and the Hessian `H`. Regularization can be applied with the optional `lasso_lambda` parameter.
   - Solve the linear system to update the projection matrix `W` and the bias `b`.

Return Results:
- After completing the specified number of iterations, return the estimated projection matrix `W`, latent variables `mu`, and bias `b`.

### Describe alternatives you've considered, if relevant

_No response_

### Additional context

In another project of mine, I have a working implementation of logistic PCA based on Mikael Brudfors's (2021) implementation under an MIT License, which I had modified to include an L1 regularization parameter per the recommendation of Lee et al. (2010). However, it is my intention, and I've already started that process, to rework the code for greater efficiency and clarity. Given that, the above-described solution will be altered as necessary while maintaining the same essential functionality. Also, I want to make sure that I don't cut Mikael out of this process, especially since to my knowledge that code was the first known implementation of this algorithm in Python. I'm seeking their involvement if they desire to be a part of this addition to the scikit-learn ecosystem if it is tentatively approved for inclusion as a new feature.

COMMENT:
@glemaitre Namely, "classical PCA" isn't meant for binary data. Similarly, "classical SVD," "classical ICA," and "classical NMF" are not meant for binary data either. Specialized logistic forms have subsequently been developed for each. 

As Song et al. (2017) discuss in their use case on genomics data:

> Compared with commonly used numerical data, binary data have some special mathematical characteristics, which should be taken into account during the data analysis. In binary measurements, '0' and '1' are abstract representations of two exclusive categories rather than numerical values 0 and 1. These two categories can also be encoded to any other two different labels, like '-1' and '1' or '-' and 'þ', without changing the meaning. Because '1' and '0' are only an abstract representation of two categories, they cannot be taken interpreted as numerical data. Furthermore, the measurement error of binary data is discrete in nature. Binary measurement error occurs when a category is assigned to the wrong label, such as when a mutated gene is misclassified as wild type. Therefore, the by default used Gaussian error assumption for continuous data in many statistical models is inappropriate for binary data analysis. Another aspect of binary data is that there can be an order in the two categories. For example, presence is often considered more important than absence. Finally, binary data can be generated from a discrete measurement process, but also a continuous measurement process [4].

>PCA is one of the most popular methods in dimension reduction with numerous applications in biology, chemistry and
many other disciplines [5]. PCA can map data points, which are in a high-dimensional space, to a low-dimensional space with minimum loss of variation. The derived low-dimensional features, which provide a parsimonious representation of the original high-dimensional data, can be used in data visualization or for further statistical analysis.

>Classical linear PCA methods are appropriate for numerical data. The direct use of linear PCA on binary data does not take into account the distinct mathematical characteristics of binary data. Although there are several parametric and nonparametric PCA alternatives for binary data in the psychometric, data analysis and machine learning literature, they are not well known to the bioinformatics community. In this article, we are going to introduce, compare and evaluate some of these different approaches.

References

Song, Yipeng, Johan A. Westerhuis, Nanne Aben, Magali Michaut, Lodewyk F. A. Wessels, Age K. Smilde. "Principal component analysis of binary genomics data." Briefings in Bioinformatics, 20.1. (2019): 317–329.

Labels: New Feature, Needs Decision - Include Feature