Linear Models from a Probabilistic Perspective: Does Logistic Regression Have an Analytical Solution?

Probabilistic Perspective

From a probabilistic perspective, mean squared error implies that we assume p(\boldsymbol{y}|\boldsymbol{x}) is a normal distribution with mean \boldsymbol{\mu}_{y|x} = \boldsymbol{x}\boldsymbol{W} + \boldsymbol{b}. Now, let’s make a stronger assumption:

Assume the joint distribution p(\boldsymbol{x}, \boldsymbol{y}) is a normal distribution.

Under this assumption, we can directly write the corresponding conditional distribution: \begin{aligned} p(\boldsymbol{y}|\boldsymbol{x}) &= \mathcal{N}(\boldsymbol{y}; \boldsymbol{\mu}_{y|x}, \boldsymbol{\Sigma}_{y|x}) \\ \boldsymbol{\mu}_{y|x} &= \boldsymbol{\mu}_y + (\boldsymbol{x} - \boldsymbol{\mu}_x)\boldsymbol{\Sigma}_{xx}^{-1}\boldsymbol{\Sigma}_{xy} \\ \boldsymbol{\Sigma}_{y|x} &= \boldsymbol{\Sigma}_{yy} - \boldsymbol{\Sigma}_{yx}\boldsymbol{\Sigma}_{xx}^{-1}\boldsymbol{\Sigma}_{xy} \end{aligned} Here \boldsymbol{\mu}_x, \boldsymbol{\mu}_y are the mean vectors of \boldsymbol{x} and \boldsymbol{y}, and \begin{pmatrix}\boldsymbol{\Sigma}_{xx} & \boldsymbol{\Sigma}_{xy} \\ \boldsymbol{\Sigma}_{yx} & \boldsymbol{\Sigma}_{yy}\end{pmatrix} is the covariance matrix of \boldsymbol{x} and \boldsymbol{y}. The form of the conditional distribution for a multivariate normal distribution can be found on Wikipedia, and its proof can be referenced on StackExchange or in relevant probability and statistics textbooks.

Comparing this with \boldsymbol{\mu}_{y|x} = \boldsymbol{x}\boldsymbol{W} + \boldsymbol{b}, we obtain: \boldsymbol{W} = \boldsymbol{\Sigma}_{xx}^{-1}\boldsymbol{\Sigma}_{xy}, \quad \boldsymbol{b} = \boldsymbol{\mu}_y - \boldsymbol{\mu}_x\boldsymbol{\Sigma}_{xx}^{-1}\boldsymbol{\Sigma}_{xy} This is, in fact, the analytical solution for linear regression.

Analysis and Reflection

Let’s summarize the process above. First, by default, linear regression only assumes the conditional distribution p(\boldsymbol{y}|\boldsymbol{x}). Through the method of least squares, we obtain the analytical solution. Then, we made a stronger assumption that the "joint distribution p(\boldsymbol{x}, \boldsymbol{y}) is a normal distribution," yet we still obtained the same analytical solution.

Why did a stronger assumption yield the same result? In fact, as seen from the loss function [eq:loss], it is quadratic with respect to \boldsymbol{y} and also quadratic with respect to \boldsymbol{x}. This means it uses at most the second-order moment information of \boldsymbol{x} and \boldsymbol{y}. Therefore, assuming the joint distribution is normal does not change the final result, as the normal distribution preserves all moment information up to the second order (mean and covariance).

Furthermore, we can imagine that the primary data statistics used by any linear model (linear regression, logistic regression, single-layer neural networks, etc.) are also moments no higher than the second order. Therefore, when dealing with linear models, we can appropriately assume a normal distribution based on the specific situation. Theoretically, this can yield an equivalent result or a sufficiently good approximation.

An Ingenious Approach

First, it is easy to verify that for the logistic regression model: \frac{p(1|\boldsymbol{x})}{p(0|\boldsymbol{x})} = \exp\left(\boldsymbol{x}\boldsymbol{w}^{\top}+b\right) \quad\Leftrightarrow\quad \ln \frac{p(1|\boldsymbol{x})}{p(0|\boldsymbol{x})} = \boldsymbol{x}\boldsymbol{w}^{\top}+b \label{eq:log} This means logistic regression is equivalent to a linear regression model with \ln \frac{p(1|\boldsymbol{x})}{p(0|\boldsymbol{x})} as the output. However, directly estimating \ln \frac{p(1|\boldsymbol{x})}{p(0|\boldsymbol{x})} is not easy. We use Bayes’ theorem: p(y|\boldsymbol{x}) = \frac{p(\boldsymbol{x}|y)p(y)}{p(\boldsymbol{x})} \quad\Leftrightarrow\quad \frac{p(1|\boldsymbol{x})}{p(0|\boldsymbol{x})} = \frac{p(\boldsymbol{x}|1)p_1}{p(\boldsymbol{x}|0)p_0} \label{eq:bys} Here p_1, p_0 are the probabilities of the positive and negative classes, which are easy to estimate. p(\boldsymbol{x}|1) and p(\boldsymbol{x}|0) are the distributions of the positive and negative samples. We now assume they follow normal distributions:

Assume p(\boldsymbol{x}|1) and p(\boldsymbol{x}|0) are normal distributions with the same covariance matrix.

Readers might find the "same covariance matrix" assumption strange; we will discuss this later. Under this assumption, let: p(\boldsymbol{x}|1) = \mathcal{N}(\boldsymbol{x}; \boldsymbol{\mu}_1, \boldsymbol{\Sigma}), \quad p(\boldsymbol{x}|0) = \mathcal{N}(\boldsymbol{x}; \boldsymbol{\mu}_0, \boldsymbol{\Sigma}) where \boldsymbol{\mu}_1, \boldsymbol{\mu}_0 are the mean vectors of the positive and negative samples, and \boldsymbol{\Sigma} can be estimated using the covariance matrix of all samples. Recalling the probability density expression for a normal distribution: \frac{1}{\sqrt{(2\pi)^n \det(\boldsymbol{\Sigma})}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})^{\top}\right\} Substituting this into Equation [eq:bys] and simplifying, we find that the quadratic terms cancel out exactly, resulting in: \ln\frac{p(1|\boldsymbol{x})}{p(0|\boldsymbol{x})} = \ln\frac{p(\boldsymbol{x}|1)p_1}{p(\boldsymbol{x}|0)p_0} = \boldsymbol{x}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0)^{\top} + \frac{1}{2}\left(\boldsymbol{\mu}_0\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_0^{\top} - \boldsymbol{\mu}_1\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_1^{\top}\right) + \ln\frac{p_1}{p_0} \label{eq:rate} Comparing this with Equation [eq:log], we get: \begin{aligned} \boldsymbol{w} &= (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0)\boldsymbol{\Sigma}^{-1} \\ b &= \frac{1}{2}\left(\boldsymbol{\mu}_0\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_0^{\top} - \boldsymbol{\mu}_1\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_1^{\top}\right) + \ln\frac{p_1}{p_0} \end{aligned} \label{eq:sol} This is an analytical solution for logistic regression. Of course, this isn’t entirely new; its logic is very consistent with "Linear Discriminant Analysis (LDA)."

Analysis and Reflection

Currently, the biggest doubt readers might have about this solution is whether the "same covariance matrix" assumption is too strong. Technically, this assumption is made so that the quadratic terms in \ln\frac{p(\boldsymbol{x}|1)}{p(\boldsymbol{x}|0)} cancel out, leaving only linear terms to obtain the analytical solution. Theoretically, is there any necessity for this assumption? In fact, we can consider that logistic regression itself (approximately) implies the assumption of the "same covariance matrix."

How to understand this? First, for the logistic regression model, Equation [eq:log] holds naturally, and Bayes’ theorem is always true. Therefore, the conclusion is that \ln\frac{p(\boldsymbol{x}|1)}{p(\boldsymbol{x}|0)} must only contain linear and constant terms. The linear regression example already told us that making a normal assumption on the data distribution of a linear model generally doesn’t lose much information. Thus, assuming p(\boldsymbol{x}|1) and p(\boldsymbol{x}|0) are normal distributions is (to some extent) reasonable. Once they are assumed to be normal, the covariance matrices must be identical for the quadratic terms to vanish.

In other words, the moment you decide to use a logistic regression model and accept the normal assumption, you have already made the assumption that "positive and negative samples have the same covariance matrix."

Multi-class Classifier

The analytical solution for logistic regression can be easily extended to the "Fully Connected + Softmax" multi-class scenario. This scenario assumes the probability of class i is: p(i|\boldsymbol{x}) = \frac{\exp\left(\boldsymbol{x}\boldsymbol{w}_i^{\top}+b_i\right)}{\sum\limits_{j=1}^k \exp\left(\boldsymbol{x}\boldsymbol{w}_j^{\top}+b_j\right)} Based on the same reasoning and assumptions, we can obtain a result similar to Equation [eq:log]: \ln \frac{p(j|\boldsymbol{x})}{p(i|\boldsymbol{x})} = \boldsymbol{x}(\boldsymbol{w}_j - \boldsymbol{w}_i)^{\top} + (b_j - b_i) And a result similar to Equation [eq:rate]: \ln\frac{p(j|\boldsymbol{x})}{p(i|\boldsymbol{x})} = \boldsymbol{x}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_j - \boldsymbol{\mu}_i)^{\top} + \frac{1}{2}\left(\boldsymbol{\mu}_i\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_i^{\top} - \boldsymbol{\mu}_j\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_j^{\top}\right) + \ln\frac{p_j}{p_i} By comparison, we find a set of solutions: \begin{aligned} \boldsymbol{w}_i &= \boldsymbol{\mu}_i\boldsymbol{\Sigma}^{-1} \\ b_i &= \ln p_i - \frac{1}{2}\boldsymbol{\mu}_i\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_i^{\top} \end{aligned}

Parameter Estimation

Finally, let’s discuss how to estimate the model parameters. We can see that \boldsymbol{w}_i, b_i are functions of p_i, \boldsymbol{\mu}_i, and \boldsymbol{\Sigma}. Thus, it essentially boils down to estimating these three.

As mentioned, p_i is simple: just use the frequency of each class. \boldsymbol{\mu}_i is also straightforward: the mean vector of each class. The difficulty lies in estimating \boldsymbol{\Sigma}. According to our previous assumption, the distribution of all data is: \tilde{p}(\boldsymbol{x}) = \sum_{i=1}^k p_i \mathcal{N}(\boldsymbol{x}; \boldsymbol{\mu}_i, \boldsymbol{\Sigma}) Integrating after multiplying both sides by \boldsymbol{x}^{\top}\boldsymbol{x}, we get: \tilde{\boldsymbol{\Sigma}} + \tilde{\boldsymbol{\mu}}^{\top} \tilde{\boldsymbol{\mu}} = \sum_{i=1}^k p_i \left(\boldsymbol{\Sigma} + \boldsymbol{\mu}_i^{\top} \boldsymbol{\mu}_i\right) = \boldsymbol{\Sigma} + \sum_{i=1}^k p_i \boldsymbol{\mu}_i^{\top} \boldsymbol{\mu}_i where \tilde{\boldsymbol{\mu}}, \tilde{\boldsymbol{\Sigma}} are the mean vector and covariance matrix of the entire dataset. Therefore, we have the estimate: \boldsymbol{\Sigma} = \tilde{\boldsymbol{\Sigma}} + \tilde{\boldsymbol{\mu}}^{\top} \tilde{\boldsymbol{\mu}} - \sum_{i=1}^k p_i \boldsymbol{\mu}_i^{\top} \boldsymbol{\mu}_i

In particular, it is recommended to perform whitening on the original data before estimation (refer to "You Might Not Need BERT-flow: A Linear Transformation Comparable to BERT-flow"), so that the mean of the entire dataset is 0 and the covariance is the identity matrix. In this case: \boldsymbol{\Sigma} = \boldsymbol{I} - \sum_{i=1}^k p_i \boldsymbol{\mu}_i^{\top} \boldsymbol{\mu}_i Theoretically, whitening before estimation is identical to direct estimation. However, for high-dimensional data, whitening is much more numerically stable, so it is recommended in practice.

Full Sample

The evaluation includes four classification tasks: Sentiment Analysis (SENTIMENT), Long Text Classification (IFLYTEK), Short News Classification (TNEWS), and E-commerce Theme Classification (SHOPPING). The general statistics are as follows:

The evaluation metric is accuracy. The results under the full training set are:

Small Sample

From the tables above, it can be seen that in terms of training set performance, the analytical solution is usually inferior to gradient descent. However, its performance on the validation and test sets is close to or even exceeds that of gradient descent. Overall, the gap between training and validation accuracy is smaller for the analytical solution, suggesting it might have better generalization capabilities. It may be more suitable for scenarios with small data volumes or where the training and validation distributions are inconsistent.

To verify this hypothesis, we kept only 1000 samples for each training set and repeated the experiments:

As seen, when training data is scarce, the gap in training accuracy narrows, but the analytical solution’s performance on the validation set comprehensively surpasses gradient descent. This further demonstrates the good generalization performance of the analytical solution in small-sample scenarios.

Comprehensive Review

These conclusions are not difficult to understand. Given that both are linear models, the analytical solution adds the assumption that "samples of each class follow a normal distribution with the same covariance matrix." When data is abundant, our estimation of each class’s distribution becomes more accurate, and the deviation of this assumption becomes more severe, making it less effective than gradient descent with adaptive training. Conversely, when data is scarce, the distribution of each class is hard to estimate anyway. In this case, the assumption acts as useful prior information, helping the model generalize "from points to surfaces," whereas gradient descent lacks such priors and suffers from insufficient generalization.

In other words, with little data, gradient descent simply memorizes a few samples and stops. The analytical solution, however, effectively "creates" more samples through additional assumptions for the model to learn from, thus learning more. With plenty of data, gradient descent "practice makes perfect" by memorizing more real samples, while the analytical solution is still "creating" samples based on its assumptions, which might be inferior to the real samples, leading to a potential drop in performance.

Is there room for improvement in the analytical solution? A direct idea is to find ways to estimate \ln \frac{p(j|\boldsymbol{x})}{p(i|\boldsymbol{x})} more precisely and then convert it into a linear regression problem to estimate parameters. There are many ways to better estimate \ln \frac{p(j|\boldsymbol{x})}{p(i|\boldsymbol{x})}, such as assuming normal distributions with different covariances or using kernel density estimation, which I leave for everyone to explore.