What's a more efficient way to calculate the max of each row in a matrix excluding its own column?

Similar idea to yours (exclude columns one by one), but with indexing:

mask = ~np.eye(cols, dtype=bool)
a[:,np.where(mask)[1]].reshape((a.shape[0], a.shape[1]-1, -1)).max(1)

Output:

array([[3, 1, 3],
       [5, 5, 2]])

You could do this using np.accumulate. Compute the forward and backward accumulations of maximums along the horizontal axis and then combine them with an offset of one:

import numpy as np

m = np.array([[1,3,1],[2,0,5]])

fmax = np.maximum.accumulate(m,axis=1)
bmax = np.maximum.accumulate(m[:,::-1],axis=1)[:,::-1]

r = np.full(m.shape,np.min(m))
r[:,:-1] = np.maximum(r[:,:-1],bmax[:,1:])
r[:,1:]  = np.maximum(r[:,1:],fmax[:,:-1])

print(r)

# [[3 1 3]
#  [5 5 2]]

This will require 3x the size of your matrix to process (although you could take that down to 2x if you want an in-place update). Adding a 3rd&4th dimension could also work using a mask but that will require columns^2 times matrix's size to process and will likely be slower.

If needed, you can apply the same technique column wise or to both dimensions (by combining row wise and column wise results).