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:
@glevv  “… over 200 known citations _between_ known papers …” It’s true that each of those referenced papers has less than 200 citations alone, but I did mention it was between them. I really wasn’t trying to shade anything there. I can provide additional sources using similar approaches. There are other papers that favorably cite Tipping (1999) in their own methods, which _alone_ have over 200 known citations. For example, Collins et al. (2001) state:

>Tipping [10] discusses a model that is very similar to our case for the Bernoulli family.

And, likewise, Schein et al. (2003) state:

>Our method adapts an algorithm, originally proposed by Tipping [15], for fitting the parameters of a closely related generative model. Like Tipping’s algorithm, it also relies on a particular convexity inequality introduced by Jaakkola and Jordan [10] for the logistic function. Our algorithm is easy to implement and guaranteed at each iteration to improve the logistic PCA log-likelihood. Thus, it provides a simple but powerful tool for dimensionality reduction of binary data.

Honestly, Lee et al. (2010) should also cite Tipping (1999) as it has citations for these other papers that both in turn do, but it does not.

References

Collins, Michael, Sanjoy Dasgupta, and Robert E. Schapire. "A generalization of principal components analysis to the exponential family." Advances in neural information processing systems 14 (2001).

Schein, Andrew I., Lawrence K. Saul, and Lyle H. Ungar. "A generalized linear model for principal component analysis of binary data." International Workshop on Artificial Intelligence and Statistics (2003): 240-247.

Labels: New Feature, Needs Decision - Include Feature