I'm working on some problems that ultimately boil down into a simple assembly of an overdetermined system of equations, $Ax=b$, where $A$ is $m \times n$ for $m \gg n$. I'm leveraging Armadillo's C++ library, which I'm linking against MKL. I computed my root mean square error on my training data by computing $E = \frac{\|Ax-b\|}{\sqrt{m}}$ and noticed that my errors and solution $x$ varied significantly per run of the test (in the second or third significant digit, sometimes some coefficients were off by a factor of 2). This troubled me, and so I ran the application with MKL_NUM_THREADS=1 forcing the number of threads down to 1, and every time I ran I got the same answer (with minor, far in the noise differences), and the error $E$ fell by several orders of magnitude.
I solve the linear system using the very simple call arma::solve(X,A,B) where $X$ contains two columns of $n$ coefficients we aim to solve for, and $B$ contains two columns of $m$ data points. The condition number is likely large (forgot to check it), but it is still within what I would guess is the "safe range" (NumPy's linalg.lstsq never complained when I tried it using that).
I'm guessing this is an error somewhere in my matrix assembly, or maybe something I'm fundamentally misunderstanding about multi-threaded MKL being used by armadillo. Is this an expected behavior or a sign that there is definitely something wrong in the code?
Thanks for any suggestions!
