# LogitModelFit crashes with large learning set

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.

• I'm surprised it works for any subset. 67 of the 784 predictors are constant (with all values being 1). Also, would you add the code that works for a smaller subset? – JimB Aug 27 '16 at 14:38
• @Jim Baldwin, I will check the constant predictors you refer to. I assume these are the fringes of the raster. Will get back on that. – Sander Aug 27 '16 at 14:45
• @Jim Baldwin, Thanks Jim, I have added the reducedData code – Sander Aug 27 '16 at 23:53

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.)

• Thanks, this improves my understanding (which is not completely there yet, but you provide valuable recommendations to digest) ... My goal is to learn the basics of machine learning. Appreciate your "black box" recommendation but it'll be for a later stage. – Sander Aug 28 '16 at 5:28
• I changed the title. – Sander Aug 28 '16 at 5:32
• Thanks, I compressed the predictor variables substantially with SVD, it now works now but still not sure if the crash is due to the number of variables or due to the rank deficiency. – Sander Sep 9 '16 at 0:45