How to speed up calculations on big symbolic matrices?

this is my first time posting something on a community of the StackExchange platform, so please feel free to correct me if I'm doing something wrong. :) Additionally you should probably know that I'm relatively new to Mathematica.

To give you some context (for tl;dr just skip below):
I'm doing some research in the field of artificial neural networks (ANN) and their properties. To be precise, I want to take a closer look at the first weight matrix to try and figure out why well-trained ANN sometimes deliver ridiculous output to seemingly proper input. To whoever might be interested in this topic: Deep Neural Networks are Easily Fooled (arxiv: 1412.1897), Intriguing properties of neural networks (arxiv: 1312.6199).
Right now I'm using one training algorithm, which is called Pseudoinverse Learning Algorithm (DOI: 10.1016/S0925-2312(03)00385-0). It's quite useful for my research as it behaves a 100% deterministically and I can build the network by just using plain matrix operations in Mathematica.

But let me get to my question:
Below you will find my code with which I generate the matrices which I want to get the Eigenvalues of. The two lines in question are the both lines defining W0. As you can see I just use the MNIST dataset included in Mathematica and prepare the data according to my needs (I only need the byte data of the pictures). The calculations run fine as long as I keep the amount of used pictures below 10 (including 10).
For 10 pictures the calculations of W0 take only half a second and the simplification of the matrix around 4 seconds. But as soon as I increase the number of the used pictures the calculation takes a considerable amount of time longer if ever it finishes - I kept my computer running for a day or two straight but the calculations didn't finish, although the sample size wasn't particularly bigger.
I already read a lot of articles regarding performance enhancements, like this one, here on StackExchange, other forums and a lot of the documentation of the involved Mathematica functions. But nothing I found and tried so far has been particularly helpful.
So my questions are, how I could speed up the calculations and if there is any explanation to the (to me seemingly odd) behavior of the long-running calculations?
I'm aware that the size of the matrix considerably grows with the size of pictures included but certainly not to an extent that a calculation with just a few more than 10 pictures fails? When I took a look around on StackExchange, some people seem to be doing calculations with matrices a lot bigger than the ones involved here. And I can't imagine that the problem only arises due to my code making use of the symbolic calculation capabilities of Mathematica.
I also know of the built-in function PseudoInverseof Mathematica which I do by hand in W0 but I want to be in control over the regularization parameter ε.

I included the remaining code in use which so far has been running fine in case the calculation of W0 finished. But I want you to be able to see what my intended use of the calculations in question is. Maybe there is a totally different (better!) approach to my problem you might want to suggest?

ClearAll["Global*"];
ClearAll[Evaluate[Context[] <> "*"]];

(* Retrieve image data from MNIST - the image download might take a while during
the first run depending on the internet connection but it
will be cached for the future *)
trainingData = ExampleData[{"MachineLearning", "MNIST"}, "TrainingData"];

stepSize = 6000;
images = trainingData[[1 ;; Length[trainingData] ;; stepSize, 1]];
getImageAsByteVector[picture_] := Flatten[ImageData[picture, "Byte"]];
imagesByteMatrix = Transpose[Map[getImageAsByteVector, images]];
Length[images]

(* After we created the data matrix over all images of MNIST we have to
add the BIAS-Neurons to the matrix *)
biasVector = Table[1., {i, 1, Length[images]}];
X0 = N[Insert[imagesByteMatrix, biasVector, 785]];

(* Create a function to dynamically insert a disturbance into the input matrix *)
X[η_] = N[ReplacePart[X0, {785, 1} -> 1. + η]];

ϵ = 10^-10;
Dimensions[X0]

minDataRange = -5.;
maxDataRange = 5.;
W0[η_] = Inverse[Transpose[X[η]].X[η] + ϵ IdentityMatrix[Length[images]]].Transpose[X[η]];//RuntimeToolsProfile
W0[η_] = Simplify[W0[η]]; // RuntimeToolsProfile
Dimensions[W0[1]]

Print["1. Square matrix W0.X0"];
WX[η_] = W0[η].X0;
WX[η_] = Simplify[WX[η]];

Print["2. Regularized square matrix σ[W0.X0]"];
workingPrecision = 16;
σWX[η_] = N[LogisticSigmoid[WX[η]], workingPrecision];

eigenvalues = N[Log[Abs[Transpose[Table[Eigenvalues[σWX[η]], {η, minDataRange, maxDataRange, 1./1.}]]]], workingPrecision]

graphs = ListPlot[eigenvalues[[1]], PlotRange -> All, Joined -> True,  DataRange -> {minDataRange, maxDataRange}]


Any help is greatly appreciated! :)
Have a nice evening and best regards!
TSwift

• I auto-formatted the symbols in your code but I think there was at least one error in the original; I changed ϵIdentityMatrix to ϵ IdentityMatrix (implicit multiplication). Please check that this is correct. – Mr.Wizard Oct 24 '16 at 7:16
• @Mr.Wizard Thank you very much! Yeah, you were right. Oddly in my code the space actually exsists, don't know what happened there. Thank you for formatting the greek letters, too! I was searching for that while I was pasting the code but couldn't find it. By the way: You seem to have a pretty close look at my code to find a missing whitespace. Any suggestions regarding my code? – TSwift Oct 24 '16 at 11:03
• Why do you need inverse on symbolic matrixes? I expect LinearSolve` would be much faster. – xslittlegrass Nov 25 '16 at 22:05