3

I am doing a comparison of some alternate linear regression techniques.

Clearly these will be bench-marked relative to OLS (Ordinary Least Squares).

But I just want a pure OLS method, no preconditioning of the data to uncover ill-conditioning in the data as you find when you use regress().

I had hoped to simply use the classic (XX)^-1XY expression? However this would necessitate using the inv() function, but in the MATLAB guide page for inv() it recommends that you use mldivide when doing least squares estimation as it is superior in terms of execution time and numerical accuracy.

However, I'm concerned as to whether it's okay to use mldivide to find the OLS estimates? As an operator it seems I can't see what the function is doing by "stepping-in" in the debugger.

Can I be assume that mldivide will produce the same answers as theoretical OLS under all conditions (including in the presence of) singular/i-ll conditioned matrices)?

If not what is the best way to compute pure OLS answers in MATLAB without any preconditioning of the data?

Joe
  • 457
  • 5
  • 18
Bazman
  • 2,058
  • 9
  • 45
  • 65

2 Answers2

3

The short answer is:

When the system A*x = b is overdetermined, both algorithms provide the same answer. When the system is underdetermined, PINV will return the solution x, that has the minimum norm (min NORM(x)). MLDIVIDE will pick the solution with least number of non-zero elements.

As for how mldivide works, MathWorks also posted a description of how the function operates.

However, you might also want to have a look at this answer for the first part of the discussion about mldivide vs. other methods when the matrix A is square.

Depending on the shape and composition of the matrix you would use either Cholesky decomposition for symmetric positive definite, LU decomposition for other square matrix or QR otherwise. Then you can can hold onto the factorization and use linsolve to essentially just do back-substitution for you.

As to whether mldivide is preferable to pinv when A is either not square (overspecified) or is square but singular, the two options will give you two of the infinitely many solutions. According to those docs, both solutions will give you exact solutions:

Both of these are exact solutions in the sense that norm(A*x-b) and norm(A*y-b)are on the order of roundoff error.

Community
  • 1
  • 1
chappjc
  • 30,359
  • 6
  • 75
  • 132
  • The pseudo code tests whether the matrix is symmetric positive definite in which case it uses cholesky and for/backward substitution, for a general square matrix is uses LU decomposition and the same for/backward substitution. Do you know which method it uses when the system is exactly determined? It's not clear cut as in the case of perfect multi-collinearity you have a de facto overdetermined system. In the pseudo code shown there is no test on the conditioning of the matrix? (Plus using "\" is the description of the \ operator is circular logic but I think I know what they mean) – Bazman Nov 01 '13 at 12:34
  • 1
    @Bazman Almost forgot about this answer, sorry. Doesn't positive definite imply exactly determined? BTW, `chol` has a second output for rank, so it may just try that first. Anyway, I also realized that I need to update my other post with at least the LU back substitution case. Thanks. The nice thing about `linsolve` is you can tell it the form of the matrix so it is efficient when inputting pre-factorized matrixes. The MathWorks should have used that instead of the circular `\` usage. – chappjc Nov 02 '13 at 15:42
1

According to the help page pinv gives a least squares solution to a system of equations, and so to solve the system Ax=b, just do x=pinv(A)*b.

David
  • 8,449
  • 1
  • 22
  • 32