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

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

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

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;

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

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