The CP tensor decomposition is used in applications such as machine learning and signal processing to discover latent low‐rank structure in multidimensional data. Computing a CP decomposition via an alternating least squares (ALS) method reduces the problem to several linear least squares problems. The standard way to solve these linear least squares subproblems is to use the normal equations, which inherit special tensor structure that can be exploited for computational efficiency. However, the normal equations are sensitive to numerical ill‐conditioning, which can compromise the results of the decomposition. In this paper, we develop versions of the CP‐ALS algorithm using the QR decomposition and the singular value decomposition, which are more numerically stable than the normal equations, to solve the linear least squares problems. Our algorithms utilize the tensor structure of the CP‐ALS subproblems efficiently, have the same complexity as the standard CP‐ALS algorithm when the input is dense and the rank is small, and are shown via examples to produce more stable results when ill‐conditioning is present. Our MATLAB implementation achieves the same running time as the standard algorithm for small ranks, and we show that the new methods can obtain lower approximation error.