R, 70 bytes

\(m,`+`=\(p)apply(p,1,\(x)convolve(c(0,x,0),(1:3)^0,,"f"))-t(!p))++m-m

Attempt This Online!

This solution is based on the built-in convolution. Convolution works on 1D vectors in R, therefore the function is applied twice (rowwise and columnwise). I came close to the @Giuseppe's answer - but unfortunately apply simplifies the output to a vector if the matrix has one single row. So the interim output has to be brought back "into shape" with -t(!p), which assures this output to be a matrix.

R, 74 73 65 bytes

\(m)mapply(\(x,y)sum(m[(n-x)^2+(p-y)^2<4]),n<-row(m),p<-col(m))-m

Attempt This Online!

-8 bytes thanks to Glory2Ukraine swapping the sapply with an mapply.

A golfed version of M--'s answer, posted with permission.

Original 79 bytes.

Explanation:

Uses the observation that an element at \$(i,j)\$ is a neighbor of an element at \$(x,y)\$ if and only if the Euclidean distance between \$(i,j)\$ and \$(x,y)\$ is less than or equal to \$\sqrt 2\$. This is then equivalent to being less than 2 since coordinates are always integers, and thus we can happily square both sides to have the squared Euclidean distance between \$(i,j)\$ and \$(x,y)\$ be less than 4. This saves 4 bytes compared to checking that \$\lvert i-x\rvert <2\$ and \$\lvert j-y\rvert <2\$.

Additionally, row(m) is guaranteed to be at least of length 1, with first element always equal to 1, so we can dispense with the ! from the usual seq(!m) by using seq(n<-row(m)) instead, which also saves another byte compared to stuffing it into the default arguments.

Wolfram Language (Mathematica), 49 48 41 bytes

-7 bytes thanks to simplification from att

ArrayFilter[#~Total~2&,#,1,Padding->0]-#&

Try it online!

Explanation

The default behavior of

ArrayFilter[f, array, r]

is to apply f to all overlapping (2r+1)×(2r+1)×… blocks of elements in array. The Padding option specifies how to handle the cases near the edge of array, where you cannot form a full block; Padding → value indicates that the "missing" elements around the edge should be filled in with value for the purposes of the calculation. For example,

ArrayFilter[f, {{1,2,3},
                {4,5,6},
                {7,8,9}}, 1, Padding → 0]

produces

In this case, I've chosen the function f to be the 2D-total of the argument, Total[#, 2] &.

But in doing so, we've summed all 9 cells in the submatrix, and we were only supposed to include the surrounding 8, not the center value itself. Thus, as a final step, subtract off the original input matrix from the result.

Dyalog APL, 30 20 bytes

{{+/,⍵}⌺(3∨1=⍴⍵)⍛-⍵}

-10 bytes thanks to @Mat_rdv who removed the need for the commute (⍨), the need to separate the operand argument (⊢), and a much more concise version of the special hanlding of 1xN matrices ((3∨1=⍴⍵)).

{⍵-⍨{+/,⍵}⌺(3 3{1⌈⍺×⍵≠1}⍴⍵)⊢⍵}­⁡​‎‎⁡⁠⁣⁤‏⁠‎⁡⁠⁤⁡‏⁠‎⁡⁠⁤⁣‏⁠‎⁡⁠⁤⁤‏⁠‎⁡⁠⁢⁣⁡‏⁠‎⁡⁠⁢⁣⁢‏⁠‎⁡⁠⁢⁣⁣‏⁠‎⁡⁠⁢⁣⁤‏‏​⁡⁠⁡‌⁢​‎⁠⁠⁠⁠‎⁡⁠⁢⁢⁡‏⁠‎⁡⁠⁢⁢⁢‏⁠‎⁡⁠⁢⁢⁣‏⁠‎⁡⁠⁢⁢⁤‏⁠‏​⁡⁠⁡‌⁣​‎‎⁡⁠⁢⁡⁣‏⁠‎⁡⁠⁢⁡⁤‏‏​⁡⁠⁡‌⁤​‎‎⁡⁠⁢⁡⁡‏⁠‎⁡⁠⁢⁡⁢‏‏​⁡⁠⁡‌⁢⁡​‎‎⁡⁠⁢⁡‏⁠‎⁡⁠⁣⁢‏⁠‎⁡⁠⁣⁣‏⁠‎⁡⁠⁢⁤⁡‏⁠‎⁡⁠⁢⁤⁢‏‏​⁡⁠⁡‌⁢⁢​‎‎⁡⁠⁢⁢‏⁠‎⁡⁠⁢⁣‏⁠‎⁡⁠⁢⁤‏⁠‎⁡⁠⁣⁡‏‏​⁡⁠⁡‌⁢⁣​‎‎⁡⁠⁣‏⁠‎⁡⁠⁤‏‏​⁡⁠⁡‌⁢⁤​‎‎⁡⁠⁢‏‏​⁡⁠⁡‌­
           (3 3{       }⍴⍵)     # ‎⁡Get the stencil shape by
                    ⍵≠1         # ‎⁢  Finding where the original input shape is not 1
                 ⍺×             # ‎⁣  multiplying by my desired shape of 3x3
               1⌈               # ‎⁤  and taking the max of (1 1) (auto-extended)
                                # ‎⁤  This is all to prevent a 1xN (or Nx1) matrix from failing on the stencil
    {    }⌺                ⊢⍵   # ‎⁢⁡Apply stencil to the input
     +/,⍵                       # ‎⁢⁢Ravel and sum
  -⍨                            # ‎⁢⁣Subtract (backwards)
 ⍵                              # ‎⁢⁤the original input matrix
💎

Created with the help of Luminespire.

I had what I thought was a more concise answer until I realized that a 3x3 stencil fails on a 1xN matrix. I'm not actually sure why; I don't understand why the normal padding rules wouldn't apply, like

3x3 stencil of 1 2 3:
┌─────┬─────┬─────┐
│0 0 0│0 0 0│0 0 0│
│0 1 2│1 2 3│2 3 0│
│0 0 0│0 0 0│0 0 0│
└─────┴─────┴─────┘

but alas, the slightly longer answer above considers this.

Python 3, 87 bytes

Probably non-competing because:

input/output is a dict that has the 2D coordinates as keys

lambda d:{(x,y):sum(v for(s,t),v in d.items()if max(abs(s-x),abs(t-y))==1)for x,y in d}

Attempt This Online!

Octave, 34 31, 28 bytes

@(x)convn(x,~~I(3),"same")-x

Try it online!

• 34 bytes: N-dimensional matrix convolution with the convenience argument "same" to have the output trimmed to input size automatically. The neighbourhood-addition-kernel matrix

  [[1 1 1]
   [1 0 1]
   [1 1 1]]
  

  is conveniently generated by inverting (~) the 3x3 diagonal matrix with only the middle element set to 1.

• 31 bytes (by @Giuseppe): z=-1:1, z|z' and cool broadcasting rules make the kernel matrix definition 3 bytes shorter.

• 28 bytes (by @Albert.Lang): convn(x,ones(3),..)-x is 2 bytes shorter than using the correct kernel, @Giuseppe remembered ~~e(3) to shorten ones another byte. I feel I(3) makes the the solution look more... complex.

R, 117 bytes

g=\(m)matrix(lapply(1:length(m),\(k){i=row(m)[k];j=col(m)[k];sum(m[abs(row(m)-i)<2&abs(col(m)-j)<2])-m[i,j]}),dim(m))

Attempt This Online!

R, 106 100 94 88 (or 70^*) bytes

Here's another method using the package {raster}:

library(raster);f=\(m)as.matrix(focal(raster(m),matrix(c(1:4,0,1:4),3,3)>0,sum,na.rm=T,pad=T))

Try this online using WebR!

_{* if we get cheeky and don't count function declaration and loading the package; for demonstration only.}

05AB1E, 18 bytes

4F¬ašøí}2Fø€ü3}OOα

Try it online or verify all test cases.

Explanation:

4F¬ašøí} # Add a border of 0s around the (implicit) input-matrix:
4F     } #  Loop 4 times:
  ¬      #   Get the first row (without popping)
         #   (using the implicit input-matrix in the first iteration)
   a     #   Convert each value to a 0 with an isLetter-check
    š    #   Prepend that row of 0s to the matrix
     øí  #   Rotate the matrix 90 degrees clockwise:
     ø   #    Zip/transpose; swapping rows/columns
      í  #    Reverse each inner row
2Fø€ü3}  # Get all overlapping 3x3 blocks:
2F    }  #  Loop 2 times:
  ø      #   Zip/transpose; swapping rows/columns
   €     #   Map over each inner list
    ü3   #    Convert it to an overlapping list of triplets
OO       # Sum each inner 3x3 block
  α      # Get the absolute difference with the values of the (implicit) input,
         # basically subtracting those values from the sums
         # (after which the result is output implicitly)

JavaScript (ES7), 84 bytes

f=(m,X,Y)=>m.map((r,y)=>r.map((v,x)=>1/X?t+=v*=(X-x)**2+(Y-y)**2<3:f(m,x,y,t=-v)|t))

Try it online!

Commented

f = (                // f is a recursive function taking:
  m,                 //   m[] = input matrix
  X, Y               //   X, Y = position of reference cell
) =>                 //   (undefined by default)
m.map((r, y) =>      // for each row r[] at index y in m[]:
  r.map((v, x) =>    //   for each value v at index x in r[]:
    1 / X ?          //     if X is defined:
      t += v *=      //       add v to t, but only if:
      (X - x) ** 2 + //         the squared Euclidean distance
      (Y - y) ** 2   //         between (x, y) and (X, Y)
      < 3            //         is less than 3
    :                //     else:
      f(             //       do a recursive call:
        m,           //         pass m[] unchanged
        x, y,        //         pass the current position
        t = -v       //         initialize t to -v
      )              //       end of recursive call
      | t            //       return t
  )                  //   end of inner map()
)                    // end of outer map()

Vyxal, 84 bits^v2, 10.5 bytes

757fbÞkḢṪƛḢṪ

Try it Online!

Bitstring:

010011011101001001000111111011000000000001011010001010100100001101100001101101010110

Convolution

Janet, 114 110 109 bytes

|(seq[[i a]:pairs $](seq[[j v]:pairs a](-(+;(seq[x :in[-1 0 1]y :in[-1 0 1]](get-in $[(+ i x)(+ j y)]0)))v)))

I always forget to write +; instead of sum…

Python 2, 105 bytes

lambda a:[[sum(sum(r[j-(j>0):j+2])for r in a[i-(i>0):i+2])-v for j,v in e(r)]for i,r in e(a)]
e=enumerate

Try it online! Link includes test cases. Explanation: Takes care not to slice before the start of the array (a[-1:2] doesn't work, although slicing past the end of the array is fine), otherwise slices the array around each value, takes the sum, and subtracts the value.

MATL, 5 bytes

3Y6Z+

Try it online!

	#implicit input, m
3Y6	#push [1 1 1; 1 0 1; 1 1 1], "neighbor cells" for convolution
Z+	#2D convolution; each element is replaced by the sum of neighbors, maintaining shape of m
	#implicit output

Charcoal, 26 bytes

ＩＥθＥι⁻ΣＥΦθ›²↔⁻κξΣΦν›²↔⁻ρμλ

Try it online! Link is to verbose version of code. Explanation:

  θ                         Input array
 Ｅ                          Map over rows
    ι                       Current row
   Ｅ                        Map over values
         θ                  Input array
        Φ                   Filtered where
           ²                Literal integer `2`
          ›                 Is greater than
               ξ            Outermost index
              κ             Inner index
            ↔⁻              Absolute difference
       Ｅ                    Map over rows
                  ν         Inner row
                 Φ          Filtered where
                    ²       Literal integer `2`
                   ›        Is greater than
                       ρ    Outer column
                        μ   Innermost column
                     ↔⁻     Absolute difference
                Σ           Take the sum
      Σ                     Take the sum
     ⁻                      Subtract
                         λ  Current value
Ｉ                           Cast to string
                            Implicitly print

Python 3, 127 bytes

lambda a:[[sum([*[0,*[[],*a,[]][i+x]][j+y:],0][0]for x in[0,1,2]for y in[0,1,2])-v for j,v in e(r)]for i,r in e(a)]
e=enumerate

Attempt This Online!

Vyxal 3, 4 bytes

Þ◌⩔∑

Vyxal It Online!

Simply return a list of all grid neighbours of each item, and sum each sublist of neighbours

Google Sheets, 48 bytes

=map(B2:D4,lambda(c,sum(offset(c,-1,-1,3,3))-c))

Replace B2:D4 with the range of cells that contain the input. To use this with variable input, add a named function with this body (44 bytes):

=map(a,lambda(c,sum(offset(c,-1,-1,3,3))-c))