QR should be more stable numerically than Cholesky on the normal equations. You don't even have to keep Q around since it will typically be large and dense (the "Q-less" QR). I presume SuiteSparseQR could be useful here?!
On Friday, July 12, 2013 10:45:14 AM UTC-7, Douglas Bates wrote: > > The good news is that I did switch to the Cholesky decomposition for > solving the weighted least squares problem at each iteration of the IRLS > algorithm and made a small change to prevent a lot of allocating/freeing > and copying in the weighting of the model matrix and these changes did save > some time. The bad news is that they didn't save a lot of time. > > I committed the cardinal sin of "optimizing" the code without first > profiling the execution. If I had done so I would have discovered that 2/3 > of the time is spent setting up the numerical representation of the model > (most of that in ModelFrame and ModelMatrix) and only 1/3 of the time in > the calculations for fitting the model. > > It is rather disappointing that even at my age there are some basic > lessons still to be learned. > > On Wednesday, July 10, 2013 5:24:49 PM UTC-5, Douglas Bates wrote: > >> On Tuesday, July 9, 2013 4:34:34 AM UTC-5, Zhuo Jia Dai wrote: >> >>> What I found useful is that in R there is a biglm package that allows >>> you to perform logistic regression on large datasets (datasets so large >>> that they can not be fully loaded into memory). I had a quick check through >>> the GLM package and it appears that it can not do this yet. Hopefully >>> someone will get around to it soon (me?) >> >> >> I have been working on a branch of the GLM repository incorporating >> methods from Dahua Lin's NumericExtensions package. In this version of the >> package I can fit a model with 53 coefficients to a binary response of >> length 2.5 million in 10.75 seconds. >> >> julia> srand(1234321) >> >> julia> df2 = DataFrame(y = float(rand(DiscreteUniform(), n)), x1 = >> rand(Normal(),n), x2 = rand(Exponential(),n), ss = >> compact(PooledDataArray(rand(DiscreteUniform(50),n)))); >> >> julia> head(df2) >> 6x4 DataFrame: >> y x1 x2 ss >> [1,] 1.0 -2.17194 0.321315 11 >> [2,] 1.0 1.1269 0.714196 16 >> [3,] 0.0 1.56257 1.52854 47 >> [4,] 1.0 0.585944 0.618314 25 >> [5,] 0.0 -1.241 0.0706758 14 >> [6,] 1.0 -1.70349 1.47446 39 >> >> >> julia> @time gm3 = glm(:(y ~ x1 + x2 + ss), df2, Bernoulli()) >> elapsed time: 10.741626201 seconds >> >> Formula: y ~ :(+(x1,x2,ss)) >> >> Coefficients: >> >> 53x4 DataFrame: >> Estimate Std.Error z value Pr(>|z|) >> [1,] 0.013152 0.00912063 1.44201 0.1493 >> [2,] 3.57473e-5 0.00126556 0.0282461 0.977466 >> [3,] -0.00313424 0.0012638 -2.48001 0.0131377 >> [4,] -0.021063 0.0127615 -1.65051 0.0988393 >> [5,] -0.00725237 0.0127732 -0.56778 0.570184 >> [6,] -0.00292751 0.0127895 -0.228899 0.818948 >> [7,] -0.0288224 0.0127543 -2.25982 0.0238322 >> [8,] -0.00808083 0.0127735 -0.632624 0.526979 >> [9,] 0.00641767 0.0127862 0.50192 0.615724 >> [10,] -0.0120857 0.012802 -0.944048 0.345145 >> [11,] 0.000269322 0.0127657 0.0210973 0.983168 >> [12,] -0.0120173 0.0127779 -0.940475 0.346974 >> [13,] -0.0133854 0.0127753 -1.04776 0.294749 >> [14,] -0.00962128 0.0127675 -0.753578 0.451102 >> [15,] -0.00964576 0.0127842 -0.754505 0.450546 >> [16,] -0.0177193 0.0127827 -1.3862 0.165686 >> [17,] -0.0158326 0.012773 -1.23954 0.215145 >> [18,] 0.0109522 0.0127804 0.856951 0.391472 >> [19,] -0.0100456 0.0127539 -0.787645 0.430905 >> [20,] -0.0176981 0.0127744 -1.38543 0.165922 >> : >> [34,] -0.0146471 0.0127598 -1.14791 0.251006 >> [35,] -0.0222798 0.0127479 -1.74772 0.0805131 >> [36,] 0.0123463 0.0127631 0.967341 0.333374 >> [37,] -0.0248778 0.0127837 -1.94606 0.0516477 >> [38,] -0.00901724 0.0127801 -0.705569 0.480456 >> [39,] -0.00358636 0.012773 -0.280777 0.778881 >> [40,] -0.00415164 0.0127966 -0.324432 0.745611 >> [41,] -0.0218191 0.0127375 -1.71298 0.0867162 >> [42,] 0.00632508 0.0127795 0.494939 0.620643 >> [43,] -0.0185398 0.0127701 -1.45181 0.146554 >> [44,] -0.00497975 0.0127506 -0.39055 0.69613 >> [45,] 0.00587257 0.0127807 0.459486 0.645885 >> [46,] -0.00246328 0.012752 -0.193168 0.846828 >> [47,] -0.0204306 0.0127699 -1.5999 0.10962 >> [48,] -0.0203154 0.0127884 -1.58858 0.112156 >> [49,] -0.00944319 0.0127884 -0.738421 0.460259 >> [50,] -0.0118946 0.0127819 -0.930581 0.35207 >> [51,] -0.00736942 0.0127909 -0.576144 0.564518 >> [52,] -0.00066287 0.0127502 -0.0519891 0.958537 >> [53,] -0.0153304 0.0127954 -1.19812 0.230872 >> >> And this isn't even the fast and space-saving version. It uses a QR >> factorization and not a Cholesky. The Cholesky should be faster. >> >> In R fitting this model with glm took 128 seconds and used up all the >> memory causing swap thrashing and wedging the system. The bigglm fit took >> 55 seconds and did not use too much memory. >> >
