def least_squares(Cui, X, Y, regularization, num_threads=<NUM_LIT:0>):
users, n_factors = X.shape<EOL>YtY = Y.T.dot(Y)<EOL>for u in range(users):<EOL><INDENT>X[u] = user_factor(Y, YtY, Cui, u, regularization, n_factors)<EOL><DEDENT>
For each user in Cui, calculate factors Xu for them
    using least squares on Y.

    Note: this is at least 10 times slower than the cython version included
    here.