3
$\begingroup$

I'm looking for a method to solve a large overdetermined system of linear equations in a least squares sense. The matrix is dense.

I'd like to use a method that works even with limited memory (we can't load the full matrix in RAM).

The matrix dimensions are something like 10,000,000 by 10,000, where I have a very large number of rows and a constant number of columns.

Improve this question
$\endgroup$
4
  • $\begingroup$ How large is the matrix in rows and columns? Although you don't have enough RAM for the full matrix, would you have enough RAM for a matrix of size n by n (where n is the number of columns)? $\endgroup$
    – Brian Borchers
    Commented Sep 25, 2014 at 14:27
  • $\begingroup$ @BrianBorchers added dimensions in question. Yes I have memory for (n,n) matix. $\endgroup$
    – mrgloom
    Commented Sep 25, 2014 at 14:37
  • $\begingroup$ What is the source problem for the matrix? Special matrix structure can be very helpful. Some dense matrices have a exploitable low-rank structure. $\endgroup$
    – Ben Thompson
    Commented Sep 25, 2014 at 22:05
  • $\begingroup$ @BenThompson I don't know structure of matrix, I'm trying to write simple general propose machine learning algorithm, where we have matrix where each row is observation and each column is variable. $\endgroup$
    – mrgloom
    Commented Sep 26, 2014 at 8:59

3 Answers 3

Reset to default
8
$\begingroup$

One option here would be to form the normal equations $A^{T}Ax=A^{T}b$ and solve them by Cholesky factorization of the resulting $n$ by $n$ matrix. This squares the condition number of the problem which could potentially be a significant problem.

Forming $B=A^{T}A$ doesn't require more than $O(n^2)$ memory, assuming that you can access the rows of $A$ one at a time. Basically,

B=zeros(n,n);

for i=1:m

B=B+A(i,:)'*A(i,:);

end

With this very large number of rows, it might be more appropriate to simply randomly sample from the rows of $A$ rather than using the entire matrix. This will of course depend on your problem data.

If ill conditioning of $A^{T}A$ is a significant problem, then you might also consider "Q-less" QR factorization methods in which you perform orthogonal transformations on $A$ to compute $R$ from the QR factorization without computing or storing $Q$ and you simultaneously perform the appropriate orthogonal transformations to the right hand side of your least squares problem.

Improve this answer
$\endgroup$
3
  • $\begingroup$ In practice, you'd probably want to load in as many rows at a time as practical and then add the outer product to $B$ using BLAS to do the work rather than working one row at a time. $\endgroup$
    – Brian Borchers
    Commented Sep 26, 2014 at 3:25
  • $\begingroup$ Anything slows down by matrix multiplication. stackoverflow.com/questions/19358984/… stackoverflow.com/questions/20983882/… Also just curious what about if I don't have space in RAM for nxn matrix? And do you know some more advanced methods? Maybe some iterative methods or randomized? I heard about randomized svd, so I think we can use pseudoinverse to solve least squares problem. $\endgroup$
    – mrgloom
    Commented Sep 29, 2014 at 9:00
  • $\begingroup$ You haven't said anything about how accurate a solution you want to have. If computing $A^{T}A$ is too slow and you're willing to settle for a less accurate solution, then you might want to use stochastic gradient descent (basically pick one equation at a time, do a gradient step with respect to that equation only and then move on to another randomly selected equation.) If you don't really need a very precise solution this might be adequate. Methods like stochastic gradient descent are widely used in machine learning because accurate solutions typically aren't needed. $\endgroup$
    – Brian Borchers
    Commented Sep 29, 2014 at 20:56
2
$\begingroup$

I think the answer to your problem is called Recursive Least Squares. Basically you treat each row one after the other, and the algorithm is based on the Woodbury matrix identity to update the inverse matrix. Enjoy ;)

Improve this answer
$\endgroup$
6
  • $\begingroup$ As far as I know I must compute initial inverse matrix at the begining anyway, or initial matrix for inverse can be (M,n) but just 1 row data and other zeros? $\endgroup$
    – mrgloom
    Commented Sep 26, 2014 at 13:36
  • $\begingroup$ Recursive least squares would be an appropriate approach to use if you need an "online" algorithm that updates the estimates as new data become available. However, the original question didn't say anything about this. $\endgroup$
    – Brian Borchers
    Commented Sep 26, 2014 at 14:23
  • $\begingroup$ @BrianBorchers : indeed, but to my opinion it remains a valid option to solve the problem. But, that is only my opinion ;) $\endgroup$
    – mellow
    Commented Sep 27, 2014 at 9:58
  • $\begingroup$ @mrgloom you can do it "online", so you can start with only one row I guess. $\endgroup$
    – mellow
    Commented Sep 27, 2014 at 9:59
  • $\begingroup$ @BrianBorchers so I can start with 1 row and reserve space for other incoming rows and fill it with zeros?(but I must define size of full matrix once at the begining? it can't be extended after?) $\endgroup$
    – mrgloom
    Commented Sep 29, 2014 at 7:49
0
$\begingroup$

For problems of the sizes that you are interested in, parallel computing seems like a good option. There are algorithms for computing QR (a key step in solving overdetermined least squares problem) in parallel computing. Mapreduce architecture is another architecture that can handle the data size that you are interested in. See this paper http://arxiv.org/pdf/1301.1071v1.pdf.

If you are limited to a single processor, you could implement the Tall skinny QR algorithm. Essentially, it tries to compute the 'R' in QR recursively and can be applied by sequentially loading parts of the matrix from disk. Details are provided in this paper http://www.netlib.org/lapack/lawnspdf/lawn204.pdf.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.