Why do we want to maximize the variance in Principal Component Analysis? I understand that in PCA, maximizing the variance is for preserving as much variability (or information) as possible during the process of reducing the dimension of the data", and I read the previous question, too: https://stackoverflow.com/questions/12395542/why-do-we-maximize-variance-during-principal-component-analysis
However, I still do not quite understand why we want to maximize the variance from the perspective of moments; for example, why don't we maximize higher-order moments, say, maximizing multiple even-order moments jointly (according to some desired weighting scheme), why just second moments?
 A: Another answer is, "we really don't care at all about maximizing variance." After all, once we get the PCs, we multiply them by 10 if we like, we rotate them, etc. For example, if the PC coefficients are very similar, like .25, .30, .27, etc, we simply re-scale the coefficients so that they are close to 1.0 and call the PC a "summate." Clearly, this destroys the variance maximization subject to unit length constraint, calling into question whether variance maximization subject to unit length constraint has any relevance.
We provide an alternative to variance maximization in our article, "Teaching Principal Components Using Correlations," recently published in Multivariate Behavioral Research, https://www.ncbi.nlm.nih.gov/pubmed/28715259 
Rather than define PCs as linear combinations that maximize variance, we summarize a (somewhat obscure) stream of literature that defines them as linear combinations that maximize average squared correlation between the linear combinations and the original variables. Then neither variance maximization nor the unit length constraint is needed. Re-scaling is allowed (encouraged even), and (non-singular) rotations are also allowed; all provide maximum average squared correlation with the original variables.
A: One answer is that maximizing variance minimizes squared error – a perhaps more immediately plausible goal.
Assume we want to reduce the dimensionality of a number of data points $\mathbf{x}_1,\cdots, \mathbf{x}_N$ to 1 by projecting onto a unit vector $\mathbf{v}$, and we want to keep the squared error small:
$$\underset{\mathbf{v}}{\text{minimize}} \, \sum_{n = 1}^N \left\|\mathbf{x}_n - (\mathbf{v}^\top \mathbf{x}_n)\mathbf{v}\right\|^2 \text{ subject to } \|\mathbf{v}\| = 1$$
This optimization problem can be turned into the equivalent problem
$$\underset{\mathbf{v}}{\text{maximize}} \,\, \mathbf{v}^\top \mathbf{C}\mathbf{v} \text{ subject to } \|\mathbf{v}\| = 1,$$
where $\mathbf{C} = 1/N\sum_{n = 1}^N \mathbf{x}_n\mathbf{x}_n^\top$. I.e., minimizing squared error is equivalent to maximizing variance along the direction of $\mathbf{v}$ (for centered data).
Another answer is that PCA is trying to fit a Gaussian model to the data (squared error and the Gaussian model are closely related). If you tried to fit another model to your data, you'd observe other moments as well (e.g., kurtosis becomes important when fitting a model via independent component analysis).
A: YAA (yet another answer).  The first distribution we encounter is the normal distribution.  In many important examples of statistics in use (linear regression, I'm looking at you) some hypothesis of normality is tacit. When statistics is taught, the normal distribution is the taught as the first and foremost example of a continuous distribution.  The connection is that the normal distribution is completely characterized by it's mean and variance, which are it's first and second moment.  So in any case of trying to understand some distribution in terms of it's moments, using any higher moment won't work sensibly if the normal distribution is among the possible outcomes. 
