tl;dr Even though this is an image classification dataset, it remains a very easy task, for which one can easily find a direct mapping from inputs to predictions.
Answer:
This is a very interesting question and thanks to the simplicity of logistic regression you can actually find out the answer.
What logistic regression does is for each image accept $784$ inputs and multiply them with weights to generate its prediction. The interesting thing is that due to the direct mapping between input and output (i.e. no hidden layer), the value of each weight corresponds to how much each one of the $784$ inputs are taken into account when computing the probability of each class. Now, by taking the weights for each class and reshaping them into $28 \times 28$ (i.e. the image resolution), we can tell what pixels are most important for the computation of each class.

Note, again, that these are the weights.
Now take a look at the above image and focus on the first two digits (i.e. zero and one). Blue weights mean that this pixel's intensity contributes a lot for that class and red values mean that it contributes negatively.
Now imagine, how does a person draw a $0$? He draws a circular shape that's empty in between. That's exactly what the weights picked up on. In fact if someone draws the middle of the image, it counts negatively as a zero. So to recognize zeros you don't need some sophisticated filters and high-level features. You can just look at the drawn pixel locations and judge according to this.
Same thing for the $1$. It always has a straight vertical line in the middle of the image. All else counts negatively.
The rest of the digits are a bit more complicated, but with little imaginations you can see the $2$, the $3$, the $7$ and the $8$. The rest of the numbers are a bit more difficult, which is what actually limits the logistic regression from reaching the high-90s.
Through this you can see that logistic regression has a very good chance of getting a lot of images right and that's why it scores so high.
The code to reproduce the above figure is a bit dated, but here you go:
import tensorflow as tf
import matplotlib.pyplot as plt
from tensorflow.examples.tutorials.mnist import input_data
# Load MNIST:
mnist = input_data.read_data_sets("MNIST_data/", one_hot=True)
# Create model
x = tf.placeholder(tf.float32, shape=(None, 784))
y = tf.placeholder(tf.float32, shape=(None, 10))
W = tf.Variable(tf.zeros((784,10)))
b = tf.Variable(tf.zeros((10)))
z = tf.matmul(x, W) + b
y_hat = tf.nn.softmax(z)
cross_entropy = tf.reduce_mean(-tf.reduce_sum(y * tf.log(y_hat), reduction_indices=[1]))
optimizer = tf.train.GradientDescentOptimizer(0.5).minimize(cross_entropy) #
correct_pred = tf.equal(tf.argmax(y_hat, 1), tf.argmax(y, 1))
accuracy = tf.reduce_mean(tf.cast(correct_pred, tf.float32))
# Train model
batch_size = 64
with tf.Session() as sess:
loss_tr, acc_tr, loss_ts, acc_ts = [], [], [], []
sess.run(tf.global_variables_initializer())
for step in range(1, 1001):
x_batch, y_batch = mnist.train.next_batch(batch_size)
sess.run(optimizer, feed_dict={x: x_batch, y: y_batch})
l_tr, a_tr = sess.run([cross_entropy, accuracy], feed_dict={x: x_batch, y: y_batch})
l_ts, a_ts = sess.run([cross_entropy, accuracy], feed_dict={x: mnist.test.images, y: mnist.test.labels})
loss_tr.append(l_tr)
acc_tr.append(a_tr)
loss_ts.append(l_ts)
acc_ts.append(a_ts)
weights = sess.run(W)
print('Test Accuracy =', sess.run(accuracy, feed_dict={x: mnist.test.images, y: mnist.test.labels}))
# Plotting:
for i in range(10):
plt.subplot(2, 5, i+1)
weight = weights[:,i].reshape([28,28])
plt.title(i)
plt.imshow(weight, cmap='RdBu') # as noted by @Eric Duminil, cmap='gray' makes the numbers stand out more
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)