By plotting the linear combinations of the principal components, it is possible to see the pixels that are being weighted in each components. Here we show the first four and the last four.

x <- mnist$train$images
plot_pc <- function(pc, title){
  a <- sweep(x, 2, colMeans(x))
  a <- rbind(a, t(pca$rotation[,pc]))
  colors = rev(RColorBrewer::brewer.pal(9, "RdBu"))

  pc_map <- apply(a, 2, function(p){
    matrix <- p[1:60000]*p[60001]
    return(matrix)
  })
  pc_map <- colSds(pc_map)
  image(1:28, 1:28, matrix(pc_map, 28, 28)[, 28:1], col = colors)
  title(title)
}

options(repr.plot.width=15, repr.plot.height=5)

layout(t(1:4))
plot_pc(1, "PC1")
plot_pc(2, "PC2")
plot_pc(3, "PC3")
plot_pc(4, "PC4")

layout(t(1:4))
plot_pc(781, "PC781")
plot_pc(782, "PC782")
plot_pc(783, "PC783")
plot_pc(784, "PC784")

By observing the figures above, we can see that the first four components cover much of the central area, where the digits are usually written, while the last four components cover very few pixels in the edges.
Before proceeding with classification models, we will apply matrix factorization with Singular Value Decomposition (SVD) to visualize the data obtained through the principal components and compare them to the full dataset.