Why isn't sum == foldl1 (+)?

Because the Data.Matrix instance for Num implements fromInteger by returning a 1x1 matrix, and + by adding elementwise which truncates the right matrix to the size of the left one (or errors if the left is larger than the right).

The sum is foldl (+) 0 which starts with a 1x1 matrix from the 0 and truncates all of the matrices in the list to that size. The foldl1 (+) does not have to make a matrix from an integer, so the result is the size of the smallest matrix in the list.

Also, sum = getSum #. foldMap Sum is the default Foldable implementation, but it is overridden by the list instance to be sum = foldl (+) 0.

The Data.Matrix implementation of + (as you show above) up to version 0.3.1.1 generates an error if the operands are not the same shape, but in version 0.3.2.0 it assumes they are the same shape without checking, and in version 0.3.3.0 it changed again (according to the comment):

-- | Perform an operation element-wise.
--   The second matrix must have at least as many rows
--   and columns as the first matrix. If it's bigger,
--   the leftover items will be ignored.
--   If it's smaller, it will cause a run-time error.
--   You may want to use 'elementwiseUnsafe' if you
--   are definitely sure that a run-time error won't
--   arise.

The implementation appears to be the same up to the current version 0.3.6.1

I think this question has been answered perfectly fine by @pat, but I will try to answer the side question:

How do I elegantly and []-safely sum up matrices?

You can't safely sum up a list of matrices, because you don't know for sure what is the size of each of the matrices in the list is, unless you have that size information at the type level. Imagine a type data Matrix (m :: Nat) (n :: Nat) e = ..., then it would be guaranteed that each of the matrices in a list have the exactly same dimensions. I am not sure if there is any work in that area for the general purpose array library, but I've seen this approach in OpenCV bindings and in mapalgebra and briefly considered using it in massiv. This approach comes with a few disadvantages, but that is not important to this question.

Another approach is to not to use a list of matrices, but instead use a 3D array, which is essentially a sequence of pages of 2D matrices of the same size. Obviously, matrix package does not have support for higher dimensions, so you can't achieve it using that library. I'll show how it can be done with massiv though. Consider we have this 3D array a:

λ> a = makeArrayR D Par (Sz3 3 4 5) $ \(i :> j :. k) -> i + j * k
λ> a
Array D Par (Sz (3 :> 4 :. 5))
  [ [ [ 0, 0, 0, 0, 0 ]
    , [ 0, 1, 2, 3, 4 ]
    , [ 0, 2, 4, 6, 8 ]
    , [ 0, 3, 6, 9, 12 ]
    ]
  , [ [ 1, 1, 1, 1, 1 ]
    , [ 1, 2, 3, 4, 5 ]
    , [ 1, 3, 5, 7, 9 ]
    , [ 1, 4, 7, 10, 13 ]
    ]
  , [ [ 2, 2, 2, 2, 2 ]
    , [ 2, 3, 4, 5, 6 ]
    , [ 2, 4, 6, 8, 10 ]
    , [ 2, 5, 8, 11, 14 ]
    ]
  ]

we can then use foldlWithin to fold across any dimension of an array, while reducing it's dimensionality:

λ> foldlWithin Dim3 (+) 0 a
Array D Par (Sz (4 :. 5))
  [ [ 3, 3, 3, 3, 3 ]
  , [ 3, 6, 9, 12, 15 ]
  , [ 3, 9, 15, 21, 27 ]
  , [ 3, 12, 21, 30, 39 ]
  ]
λ> foldlWithin Dim1 (+) 0 a
Array D Par (Sz (3 :. 4))
  [ [ 0, 10, 20, 30 ]
  , [ 5, 15, 25, 35 ]
  , [ 10, 20, 30, 40 ]
  ]

Note that there is no chance of an exception to be raised or some other unpredictable behavior like you've experienced at any point here (putting async exceptions aside, of course). Alternatively it is possible to take slices of the array, and add those individually, but that approach is a bit more involved if we want to avoid exceptions:

λ> a !?> 0
Array D Par (Sz (4 :. 5))
  [ [ 0, 0, 0, 0, 0 ]
  , [ 0, 1, 2, 3, 4 ]
  , [ 0, 2, 4, 6, 8 ]
  , [ 0, 3, 6, 9, 12 ]
  ]
λ> import Data.Maybe
λ> fromMaybe empty $ a !?> 0 >>= \acc0 -> foldlM (\acc pageIx -> (acc +) <$> (a !?> pageIx)) acc0 (1 ... 2)
Array D Par (Sz (4 :. 5))
  [ [ 3, 3, 3, 3, 3 ]
  , [ 3, 6, 9, 12, 15 ]
  , [ 3, 9, 15, 21, 27 ]
  , [ 3, 12, 21, 30, 39 ]
  ]

Above will produce Nothing if a doesn't have a page with index 0 through 2. On the other hand, if you are ok with runtime exceptions, then this will be a bit clearer, but not as safe:

λ> foldlS (\acc pageIx -> acc + (a !> pageIx)) (a !> 0) (1 ... 2)
Array D Par (Sz (4 :. 5))
  [ [ 3, 3, 3, 3, 3 ]
  , [ 3, 6, 9, 12, 15 ]
  , [ 3, 9, 15, 21, 27 ]
  , [ 3, 12, 21, 30, 39 ]
  ]

Or using lists, as it was stated in the actual question:

λ> Prelude.foldl1 (+) $ Prelude.fmap (a !>) [0 .. 2]
Array D Par (Sz (4 :. 5))
  [ [ 3, 3, 3, 3, 3 ]
  , [ 3, 6, 9, 12, 15 ]
  , [ 3, 9, 15, 21, 27 ]
  , [ 3, 12, 21, 30, 39 ]
  ]

A side note. I noticed the usage of the SECIALIZE pragma, which leads me to believe that you actually care about performance. A small, but important fact on that: matrix package uses boxed vectors underneath for Matrix datatype, which will always give you terrible performance and no pragma can help you with that.