PCA gives me loadings with different signs. I understand that I may simply revert them (as explained in this thread), but I need an explanation (my boss insists) of why Matlab implementation gives me these signs, but not reverted.
UPDATE I found empirically, that in Matlab the loading with highest absolute value is always positive.
Given any solution to PCA, a sign-flipped version of it is an equally valid solution. A numerical solver breaks this symmetry by finding one of these equally valid solutions. Implementation details and initial conditions determine which solution the solver will produce.
One way to think about PCA is that it maximizes the sum of the variance of the data projected onto the weight vectors, subject to the constraint that the weight vectors are orthonormal. Say the data set contains $n$ points in a $d$ dimensional space with mean $\mu$. We seek a set of orthonormal vectors $\{v_1, ..., v_p\}$ that solves the following optimization problem:
$$\max_{v_i, ..., v_p} \quad \sum_{i=1}^p \frac{1}{n} \sum_{j=1}^n \left [ (x_j - \mu)^T v_i \right ]^2 \quad \quad s.t. \quad \begin{array}{ll} \|v_i\| = 1 & \forall i \quad \\ v_i^T v_j = 0 & \forall i \ne j \\ \end{array} $$
Say we have a set of weight vectors that solves this problem. You can see that flipping their signs gives the exact same value for the objective function, and the constraints remain satisfied. Hence, the sign-flipped solution is equally valid. There are various other ways of thinking about PCA and writing it as an optimization problem. Some of these ways sound conceptually different, but they all yield the same set of solutions, and the same reasoning holds for all of them.
The reason a particular function implementing PCA returns any one of these solutions over the others comes down to implementation details and initial conditions. Say we're trying to solve the above problem using a standard optimization solver. It starts from some (possibly random) initial set of parameters, then iteratively updates them to increase the value of the objective function, while respecting the constraints. Imagine the objective function as a hilly landscape, where each location corresponds to a particular choice of parameters and the height at each location is the value of the objective function for those parameters. The constraints define particular regions the solver is allowed to go. Each solution to the problem is the highest allowed location on some surrounding hill. There are multiple hills, and the solutions all have the same height (i.e. are equally good). The solver starts from some initial location in this landscape and generally tries to move uphill, eventualy stopping when it can't make any further uphill progress. So, the solution it finally attains is determined by the hill it starts on and how it chooses to step around the landscape.
Of course, one wouldn't typically solve PCA this way because there are more specialized, computationally efficient ways to do it. For example, one popular method is to obtain the weights as eigenvectors of the covariance matrix. But the eigenvalue solver is itself an iterative algorithm, and is subject to the same kinds of issues.
An eigenvector remains an eigenvector when you change its sign. Thus, an eigenvector does not have a "sign". As loadings as just eigenvectors multiplied by the square root of the associated eigenvalue, they have no "sign" either.
So now, why MAtlab gives you a specific sign I guess you would have to see how the decomposition is made behind. As far as I have seen on different topics (there or there) the sign is chosen arbitrarily.
Maybe this can help you. Otherwise, there has been some attempt to define what a "right" sign would be here.
EDIT :
Here is a website listing the methods chosen to solve the problem depending on the type of input (type of data in the matrix). I did not want to put it as I cannot find the sources the author of the page uses. But as far as I have seen, it seems to be still true:
