Re: least square solution of b = A * x with a very huge A matrix

[c.]

On 31 May 2013, at 00:25, andrea console <address@hidden> wrote:

> Thank you all.
> In this example the matrices size is:
> A: 1072452 X  323136 (sparse)
> b:  1072452  X 1
> c:  323136 X 1
what is c?
> This is the file (67MB): 
> https://docs.google.com/file/d/0B4nOH7xkcKiabkdSNktLT0xhYjg/edit?usp=sharing
> 
> I didn't save it as blocks because, before computation, I have to remove some 
> rows (and corresponding elements in b array), so the block structure is not 
> preserved
> 


The matrix is very large but very sparse and "very rectangular". 
As the original block structure is lost, the best I can suggest at this point is

AA = A' * A;
bb = A' * b;

%% then either
x = AA \ bb; %% this will take veeeeery long
%% or
x = pcg (AA, bb, [], []); 
%% this will be inaccurate but you can refine the solution by an additional step
x = pcg (AA, bb, [], [], x); %% repeat ad libitum

Alternatively you may  want to use a singular value decomposition, unfortunately
I just found a small bug in svds.m that creates a full size vector from the 
sparse matrix input
at an intermediate step of the algorithm, which I fixed with this changeset:

http://hg.savannah.gnu.org/hgweb/octave/rev/bd5e1b3a1ef9

so you'll also need to apply the same change to svds.m before you can use it:
 
--- a/scripts/sparse/svds.m 
+++ b/scripts/sparse/svds.m
@@ -136,8 +136,9 @@

      error ("svds: SIGMA must be a positive real value or the string 'L'");
    endif

-  [m, n] = size (A);
-  max_a = max (abs (A(:)));
+  [m, n] = size (A);  
+  max_a = max (abs (nonzeros (A)));
+
    if (max_a == 0)
     s = zeros (k, 1);  # special case of zero matrix
    else

c.