How to solve a set of linear systems

I’m having a problem understanding the working mechanism of tensorflow’s function: tf.linalg.solve. I want to solve a set of linear systems (AX = Y), where the linear coefficients (A) were shared but there are multiple batches of Y, which are different. Using numpy, I can simply do it via:

np.random.seed(0)
mtx = np.random.uniform(size= (1,4,4))
vec = np.random.uniform(size= (100,4,1))
solution = np.linalg.solve(mtx,vec)
print(abs(np.matmul(mtx,solution) - vec).max())
# 5.551115123125783e-16

which gives me a quite consistent solution. But when I switch to tensorflow, it gives me the results:

mtx = tf.random.uniform(shape = (1,4,4))
vec = tf.random.uniform(shape = (100,4,1))
solution = tf.linalg.solve(mtx,vec)
print(tf.math.reduce_max(abs(tf.matmul(mtx,solution) - vec))) 
# tf.Tensor(1.3136615, shape=(), dtype=float32)

According to the document, I assume the solution should be solved according to the corresponding vec. But it does not seem to give me the expected results in tensorflow. Since I’m a new user, I could have messed up something. It would be appreciated if any information could be offered.