Efficient discrete Laplacian of a matrix

ListCorrelate is a very general function to handle operations of this type. The Laplacian is a correlation with a specific kernel.

ker = {{0, 1, 0}, {1, -4, 1}, {0, 1, 0}};
ListCorrelate[ker, Z]

which gives the same answer as your laplacian function. As for speed, this appears to be faster than laplacian for small sizes, slower for medium sizes, and faster for larger sizes. Consider the comparison:

n = 2^14;
mat = RandomReal[10, {n, n}];
AbsoluteTiming[laplacian[mat]][[1]]
ker = {{0, 1, 0}, {1, -4, 1}, {0, 1, 0}};
AbsoluteTiming[ListCorrelate[ker, mat]][[1]]
(*
  81.0357
  35.7189 
*)

Your compiling fails because 3;; can't be compiled, see here for more information. The easiest fix is to modify it to 3 ;; -1, but this leads to no improvement for speed. To write a faster laplacian, we can make use of the experience obtained in this post:

cLa = Hold@Compile[{{z, _Real, 2}}, 
     Module[{d1, d2}, {d1, d2} = Dimensions@z; 
      Table[z[[i + 1, j]] + z[[i, j + 1]] + z[[i - 1, j]] + z[[i, j - 1]] - 
        4 z[[i, j]], {i, 2, d1 - 1}, {j, 2, d2 - 1}]], CompilationTarget -> C, 
     RuntimeOptions -> "Speed"] /. Part -> Compile`GetElement // ReleaseHold;

mat = RandomReal[10, {1000, 1000}];
Array[laplacian[mat] &, 10]; // AbsoluteTiming
(* {0.447978, Null} *)
Array[cLa[mat] &, 10]; // AbsoluteTiming
(* {0.084375, Null} *)

You can get very good speed and additional flexibility by using the built-in functionality described in "The Numerical Method of Lines". For example, here I assume that the dimensions of the matrix you call Z are known a priori, so that they can be hard-coded. That allows me to define a matrix d2 that performs the Laplacian on the flattened version of Z:

d2 = 
  SparseArray@
   N@Sum[NDSolve`FiniteDifferenceDerivative[i, {#, #} &[Range[1000]], 
       "DifferenceOrder" -> 2]["DifferentiationMatrix"], {i, {{2, 0}, {0, 2}}}];

lap2[values_] := Partition[d2.Flatten[values], Length[values]]

Z = RandomReal[10, {1000, 1000}];
Array[lap2[Z] &, 10]; // AbsoluteTiming

(* ==> {0.079342, Null} *)

In the function FiniteDifferenceDerivative, I could now also specify higher differentiation orders than 2. The generalization to three dimensions would also be straightforward, by changing the sum to run over {{2,0,0},{0,2,0},{0,0,2}} instead of {{2, 0}, {0, 2}}.

Another option is to add ,PeriodicInterpolation->True after the "DifferenceOrder" option. So this lets you do lots of variations on the Laplacian very efficiently.