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.
pca()
several times on the same data matrix, one gets the same signs (at least usually). I suppose that it has little to do with Matlab, but depends on some specific initialization choices in ARPACK. $\endgroup$ – amoeba Feb 1 '17 at 19:35