calling all "R Statistical Package" gurus

[Ether]

Given m-by-n matrix [A], (m>n), and m-by-1 column vector [b] for the overdetermined system of linear equations

[A][x] = [b],

What is the recommended way to use R to find the n-by-1 column vector [x] which minimizes the L1 norm of [b]-[A][x] ?

aka minx |[b]-[A][x]|₁, aka Least Absolute Deviation (LAD) ?

[faust1706]

This paper might be of use to you.

I have personally used rgenoud, but my code for it is sitting on the desktop in my lab at university and wouldn't be much use to you because it maximize the cost function described in this paper.

[GeeTwo]

Whew. I somehow managed to get an MS in physics without taking a single linear algebra course (I don't know how, and I don't advise it!). When I've had this sort of problem since entering the workforce, I just trusted MATLAB to figure out the best method.

[Ether]

Originally Posted by GeeTwo

Whew. I somehow managed to get an MS in physics without taking a single linear algebra course (I don't know how, and I don't advise it!). When I've had this sort of problem since entering the workforce, I just trusted MATLAB to figure out the best method.

I already have an Octave (Matlab) script that does this. This thread is about the recommended way to do it using R.

[Ether]

I found this statement:

Package quantreg contains variations of simplex and of interior point routines ( nlrq(), crq()). It provides an interface to L1 regression in the R code of function rq(). [SPLP, LP, IPM]

... at this webpage:

http://cran.r-project.org/web/views/Optimization.html



I also found this statement:

Koenker's quantile regression package quantreg contains L1 (aka LAD, least absolute deviations)-regression as a special case

... at this webpage:

http://cran.cnr.berkeley.edu/web/views/Robust.html


Attached is a zip file containing an [A] matrix and [b] column vector in CSV format. Would some kind R guru please show how to use the above method to do L1 regression on the attached data? (see original post in this thread)

[Ether]

Originally Posted by Ether

Given m-by-n matrix [A], (m>n), and m-by-1 column vector [b] for the overdetermined system of linear equations

[A][x] = [b],

What is the recommended way to use R to find the n-by-1 column vector [x] which minimizes the L1 norm of [b]-[A][x] ?

I found an elegant one-line-of-source-code fast-executing solution.

I knew R should be able to do this.

In terms of code simplicity (using tested and stable publicly available supported free libraries) I would rank R at the top of the list (for this particular problem), tied for first place with Python. In second place is AMPL. Matlab/Octave is a distant third.

PM me if you're interested.


PS: I haven't looked into SciLab yet. How about LabVIEW?? Any takers?

[Ether]

Originally Posted by Ether

I would rank R at the top of the list (for this particular problem)

I should add the following qualifier: the above is true for small datasets (like analyzing one event). For larger datasets (like doing world OPR for an entire season's worth of events), not so much.

Originally Posted by Ether

Matlab/Octave is a distant third.

Matlab/Octave has good support for sparse matrices to effectively process larger datasets.

[Rachel Lim]

Originally Posted by Ether

In terms of code simplicity (using tested and stable publicly available supported free libraries) I would rank R at the top of the list (for this particular problem), tied for first place with Python. In second place is AMPL. Matlab/Octave is a distant third.

Is the difference just with how easy it is to write the code / input the data, or do they do the calculations differently? If the latter is true, would different CAS programs fall into a few different groups, or would they all do the calculations differently?

[Ether]

R uses the rq() function in the quantreg package which does quantile regression. It converts the problem to SDP and solves it. If you set the tau parameter of rq() to 50%, it will do L1 regression.

Python has a convex optimization package that contains a function for L1 regression. You don't have to know how to setup the problem, just pass it your A & b data.

AMPL is a modeling language which allows you to state the problem as a constrained minimization, and then AMPL expands it into a linear program format that can be passed to your linear solver of choice.

Octave has glpk, the GNU Linear Programming kit. You have to manually code the problem as a linear program and then pass it to glpk() for solution.

Matlab has the linprog() function which does the same thing as Octave's glpk(). It also has the CVX package, which allows you to use Matlab as a modeling language, but it is not intended for large problems.