Is PCA a continuous function of the data? Suppose that my data are such that a PCA gives a unique solution for the first principal component up to scaling (e.g. my data do not all lie on a circle, or some such weirdness).
Is it the case that the first principal component is continuous as a function of the data $X$?  I.e. if we move one of the $x \in X$ by a small amount $\epsilon$, then the direction of the first principal component changes by a small amount?
If the function is discontinuous, what are the constraints within which it is a continuous function?
EDIT:
To try and clarify, with some assumptions barring numerical edge cases (some of which are enumerated here), we can think of PCA$_1$ as a function that takes in some data, $X$ of $n$ observations in $R^m$ (i.e. $X \in R^{n\times m}$), and gives back the direction of the first principal component PCA$_1(X)$.  The actual vector is unique up to a scalar multiple, so to be formal we can reduce it to mere direction by, say, defining PCA$_1$ to hand back a vector where the first variable is positive (and if it is zero, then where the next variable is positive, and so on) of unit length.  Consider some alternate matrix $X' \in R^{n\times m}$.  Is it the case that for every $\delta > 0$ there exists an $\epsilon > 0$ so that if $||X-X'||_F < \epsilon$ then $||$PCA$_1(X) - $PCA$_1(X')||_F < \delta$?
 A: I don't think it is. Imagine the data that is almost spherical with covariance matrix $[[1+\epsilon,0],[0,1]]$ and that after a change of one vector becomes $[[1,0],[0,1+\epsilon]]$. Then the first component will flip by 90 degrees.
A: Definitions
According to the Wikipedia article, a continuous function is a function where small changes in the inputs lead to arbitrarily small changes in the output. A function can be discontinuous but within certain limits be continuous.
PCA is the decomposition of data into a basis set of orthogonal basis vectors (basically the individual building blocks of the dataset) that explain the maximal amount of variance present in the dataset. It is based on the assumption that dataset $D$ can be represented by a sequence of simultaneous equations and the dataset is the product of scale matrix $S$ and the transpose of the simultaneous equations/basis vectors $L^T$. $$ D = SL^T$$
These basis vectors explain different proportions of the variance in the dataset, and by convention we rank these in order of descending contribution as evaluated by the eigenvalue. Numerous algorithms exist to calculate PCA, some (such as NIPALs) explicitly calculate in ranked order by iteratively calculating PCs and imposing a maximisation constraint. The mathematical identity of a PC is not its rank, but its shape/direction with respect to the original data.
So I see two parts that need considered to answer the question; ‘Is the implementation of PCA a continuous function?’ and ‘Is PCA itself a continuous function?’
I’ll draw on an example frommy distance past, looking at factors influencing nutritional quality of milk. In this example our dataset is influenced by 4 independent factors, $BREED$, $FEED$,$LACTATION$,$SEASON$ and we want to explore the impact of perturbation on any one these factors (let’s choose $BREED$). None of our variables are direct measures of these underlying factors, in fact it is a nutritional profile of the milk. Our PCA returns 4 large eigenvalue PCs after which the eigenvalues drop dramatically. We can determine that each of the PCs are largely accounted for by the different factors because we see strong correlations with specific PCs and the values of specific underlying factors.
Is PCA a continuous function?
Calculation of PCA is based on OLS of the data by itself. We calculate the sum of the variance (sum of squares, this is done by inner product of the data matrix $D$). For the $i^{th}$ PC:
$$ l_i^T = argmin(|D_{i-1}^TD – s_i^+D{i-1}|) $$
We then identify what combination of weights $s_i$ applied to the data allows us to maximise our estimation of the sum of squares (minimise the residual). This means that for the calculation of PCA to be discontinuous we need to identify a discontinuity in the sum function, the square function or the combination of the two. Lets say we find that $BREED$ is the largest contributor to the variance observed in the dataset, with $FEED$ a close second. First we choose to limit the variance perturbation to less than the difference in variances (square of the eigenvalue) for the two factors ($\Sigma \epsilon^2 <\Sigma BREED^2 - \Sigma FEED^2$), then the perturbation is sufficiently small to prevent $FEED$ usurping $BREED$ as the first ranked PC.
If we focus on  the first PC in our original dataset then perturb our dataset by a small value $\epsilon$ the this becomes
$$ l_1^T = argmin(|(D+\epsilon)^T(D+\epsilon) – s_1^+D+\epsilon|) $$
This means that both $l_1^T$ and $s_1$ are perturbed by an amount that satisfies the argmin condition. The impact of $\epsilon$ depends on its covariance with the original dataset.
If it is perfectly positively correlated, then all of the variable magnitudes in $\epsilon$ fit the original $l_1^T$, meaning its shape is unchanged. However, the variance of the perturbation is now added to the original variance of this PC, increasing its eigenvalue. Since the definition of PCs is orthogonality, if it is perfectly correlated with $l_1^T$ then it will provide 0 variance to the other PCs. In this scenario the impact on the shape of any $L_T$ is zero, but the eigenvalues are modified for the correlated PC. The sum of squares is a continuous function.
If it is perfectly negatively correlated, then all of the variable magnitudes in $\epsilon$ oppose the original $l_1^T$, meaning its shape is unchanged because $\epsilon$ is much smaller than the data for continuity testing. However, the variance of the perturbation is now removed from the original variance of this PC, decreasing its eigenvalue. Since the definition of PCs is orthogonality, if it is perfectly correlated with $l_1^T$ then it will provide 0 variance to the other PCs. In this scenario the impact on the shape of any $L_T$ is zero, but the eigenvalues are modified for the correlated PC.
If it is orthogonal then none of the variable magnitudes fit the original $l_1^t$ and it reduces the fit quality wrt the sum of squares. This means that a new $l_1^t$ will be calculated, but $\epsilon$ will only disturb the sum of squares by $\epsilon^2$, so if $\epsilon$ tends to zero, its square tends to zero quicker. This means for small perturbations by orthogonal variables we would see it have an impact proportional to its square. The square function is a continuous function.
If the correlation is between 0 and 1 then we see a balance between strengthening the existing PC (increasing its eigenvalue by the magnitude of $\epsilon$) and introducing noise (changing its shape by the square of the correlation moment between $\epsilon$ and $l_1^T$). Correlation is a continuous function between its limits of -1 and 1.
I would suggest that the above suggests that small continuous changes in the data induce small continuous changes in the vectors calculated as long as the above assumption about the magnitude of $\epsilon$ is met. However, it is important to realise that the perturbation is not just introducing changes in variance, it is also introducing changes in covariance, and so where exactly the perturbation becomes evident is dependent on the level of covariance. It may be a bit messy but there is no discontinuity in this behaviour.
Next we consider where this may break down
Is the implementation of PCA a continuous function?
If we follow the conventional process and specifically select the highest ranked eigenvalue, where does the continuity break down? The breakdown will occur when a lower ranked PC switches rank (this is the almost spherical situation described in the question and expanded on by @sega_sai) – leading to PC1 being returned as a completely different underlying factor, in our example this would likely be $FEED$ being returned at the expense of $BREED$.
When does this occur?
If $\epsilon$ is perfectly correlated with $FEED$ then their covariance will be $\Sigma \epsilon^2 + \Sigma FEED^2$ This means that if $\Sigma \epsilon^2 =\Sigma BREED^2 - \Sigma FEED^2+1/\infty$ then ranks will switch.
More generally if $\Sigma FEED^2 +\Sigma FEED*\epsilon>\Sigma BREED^2 +\Sigma BREED*\epsilon$ then we will observe a rank switch – i.e. if the sum of the variance of $FEED$ and its covariance with $\epsilon$ exceeds the sum of variance of $BREED$ and its covariance with $\epsilon$ then rank will switch. This can be achieved by increasing correlation with $FEED$ at the expense of $BREED$ or by having a weak correlation and a massive variance for $\epsilon$.
However, an important feature of the PCA equation is that rank or order do not matter, it is by convention (and for some algorithms practical implementation) that the vectors are ordered by descending eigenvalue.
Instead of arbitrarily demanding that rank act as an identifier of PCs, we can choose to instead define them by their shape/direction. This means we perform a correlation analysis between the PCs for the original and perturbed models. We use the correlation to find the PCs which are most similar and so it will match the PCs no matter what rank. This would avoid the discontinuity associated with rank based selection. The correlation will vary as a continuous function, dependent on the nature of the interaction between $\epsilon$ and the original data.
