LogitModelFit crashes with large learning set

Question

I am using LogitModelFit to implement a classification system for handwritten digits, based on the MNIST ExampleData set. The problem is that the full learning set - A 60,000 x 784 FP matrix - crashes LogitModelFit. A smaller set works, but I was wondering how to avoid the crash as I (right or wrong) assumed it should be capable handling these size learning sets.

I know this can be done with Classify, but it's a bit too black box for me and I am not sure if does cope with similar size trainingsets.

This is what I do:

Import the data

{X, Y} = (ExampleData[{"MachineLearning", "MNIST"}, "Data"] /. Rule -> List)\[Transpose];

Extract imagedata, vectorize the raster and define vector size

θsize = AbsoluteTiming[Dimensions[matrixX =Flatten[Map[ImageData, X], {3, 2}]\[Transpose]]][[2,2]];
θ= Table[Symbol["θ" <> ToString[i]], {i,θsize}];

Separate training and test data (X) and predicates (Y)

trainingX = matrixX[[;; 60000]];
trainingY = Y[[;; 60000]];
testX = matrixX[[60001 ;;]];
testY = Y[[60001 ;;]];
Clear[X, Y, matrixX]

Identify ranges for digits 0..9

individualTrainingSets = MovingMap[List[#[[1]] + 1, #[[2]]] &, Prepend[Accumulate[Length /@ Split[trainingY]], 0], 2]

(* {{1, 5923}, {5924, 12665}, {12666, 18623}, {18624, 24754}, {24755, 30596}, {30597, 36017}, {36018, 41935}, {41936, 48200}, {48201, 54051}, {54052, 60000}}

Yranges = Table[Symbol["trainingY" <> ToString[i]], {i, 0, 9}];
Y = SparseArray[i_ /; #1 <= i <= #2 -> 1, {60000}] & @@@individualTrainingSets;

Construct trainingset for "0"

data = Transpose[trainingX\[Transpose]~Join~{Y[[1]]}];

This work for smaller subsets, EDIT

reducedData = Take[data, {1, 60000, 100}];

rlm = LogitModelFit[reducedData,θ,θ]; // AbsoluteTiming

END EDIT

crashes for the full set (without error message)

rlm = LogitModelFit[data,θ,θ]; // AbsoluteTiming

I would like to understand why this function cannot cope with the size of the training set and - possibly - how to resolve this.

Answer 1

First, this example data is not associated with LogitModelFit so I don't think you should have the title of the question suggest that one of LogitModelFit's examples doesn't work. Second, I guess I have low expectations because I can't conceive of using any generalized linear model with 784 parameters in a single model. (Looking at multiple subsets of 784 predictor variables, yes, but all at once, no.)

At most you can have is 713 of the predictor variables because that is the rank of the design matrix:

x = data;
(* Replace the dependent variable with a constant for the intercept
 although it doesn't really matter for this example *)
x[[All, 785]] = Table[1, {i, 60000}];
designMatrix = Transpose[x].x;
MatrixRank[designMatrix]
(* 713 *)

Update

One of the difficulties with using the example dataset with logistic regression is that 231 of the 784 variables have less than 100 unique values out of 60,000 observations and the maximum number of unique values for a variable is 256.

I take it that the underlying question isn't really about the example dataset but rather how large a data set can LogisticModelFit fit. It can certainly handle 60,000 observations and 400 predictor variables. Here I've kept the 400 predictor variables with at least 255 unique values:

keep = {130, 131, 148, 150, 151, 152, 154, 155, 156, 157, 158, 159, 
   160, 161, 162, 163, 175, 176, 177, 178, 179, 180, 181, 182, 183, 
   184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 202, 203, 204, 
   205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 
   218, 219, 220, 221, 230, 231, 232, 233, 234, 235, 236, 237, 238, 
   239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 257, 
   258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 
   271, 272, 273, 274, 275, 276, 277, 278, 285, 286, 287, 288, 289, 
   290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 
   303, 304, 305, 306, 313, 314, 315, 316, 317, 318, 319, 320, 321, 
   322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 
   340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 
   353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 368, 369, 370, 
   371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 
   384, 385, 386, 387, 388, 389, 390, 396, 397, 398, 399, 400, 401, 
   402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 
   415, 416, 417, 418, 424, 425, 426, 427, 428, 429, 430, 431, 432, 
   433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 
   446, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 
   464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 480, 481, 
   482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 
   495, 496, 497, 498, 499, 500, 501, 508, 509, 510, 511, 512, 513, 
   514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 
   527, 528, 529, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 
   546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 565, 566, 
   567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 
   580, 581, 582, 583, 584, 594, 595, 596, 597, 598, 599, 600, 601, 
   602, 603, 604, 605, 606, 607, 608, 609, 610, 611, 621, 622, 623, 
   624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 
   637, 638, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 
   661, 662, 663, 664, 665, 666, 679, 680, 681, 682, 690, 785};
data2 = data[[All, keep]];
θ = Table[ToExpression["θ" <> ToString[keep[[i]]]], {i,  Length[keep] - 1}];
Dimensions[θ]
rlm = LogitModelFit[data2, θ, θ]; // AbsoluteTiming
(* {456.929, Null} *) 

For most software packages the fitting of such models has the number of predictors being a far more limiting factor than the number of observations. But when one is mainly or solely interested in "predicting" rather than "explaining" with zillions of predictor variables, you're probably better off going the black box route (but checked with appropriate cross-validation, out-of-bag error estimates, etc.)