I think you can use 2D convolution to find the places where the small array b can go in the large array a. If you use scipy.signal.convolve2d with mode='valid' you only get locations where the small array 'fits'. I think using the abs of the arays gets around positive and negative values (in either array) canceling, but I haven't tested any of this very rigorously.

Here's what I did, using @CypherX's fill_a_with_b function for the fill step:

import numpy as np
import scipy.signal

# Your input data.
a = np.array([[0, 0, 0, 9], 
              [0, 0, 0, 7],
              [0, 0, 0, 2],
              [2, 3, 1, 5]])
b = np.ones((2, 2)) * 3

# Find places where b can go.
allowed = scipy.signal.convolve2d(np.abs(a), np.abs(b), mode='valid')

# Get these locations as (row, col) pairs.
coords = np.stack(np.where(allowed==0)).T

# Choose one of the locations at random.
choice = coords[np.random.randint(coords.shape[0])]

# Use @CypherX's 'fill' function.
def fill_a_with_b(a, b, pos=[0, 0]):
    aa = a.copy()
    aa[slice(pos[0], pos[0] + b.shape[0]), 
       slice(pos[1], pos[1] + b.shape[1])] = b.copy()
    return aa

# Do the fill thing.
fill_a_with_b(a, b, choice)

This results in (for example)...

array([[0, 0, 0, 9],
       [0, 3, 3, 7],
       [0, 3, 3, 2],
       [2, 3, 1, 5]])

I will try giving you an example. But it will be based on some assumptions:

1. You know that the 0's are spanned across a contiguous rectangular block.
2. There are no other zeros in a.

If you would like to fill in a non-contiguous block of zeroes or, there are zeros on the columns/rows where you have some other non-zero values, you would have to think of a more sophisticated solution.

Solution: Random Insertion of array b into array a where a==0

Assumption: we know that the places where a is zero are a contiguous set of positions with a rectangular shape.

Imports

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
%config InlineBackend.figure_format = 'svg' # 'svg', 'retina'
plt.style.use('seaborn-white')

Make Data

# Make a
shape = (5,5)
a = np.zeros(shape)
a[:,-1] = np.arange(shape[0]) + 10
a[-1,:] = np.arange(shape[1]) + 10
# Make b
b = np.ones((2,2))*2

Preprocess

Here we determine possible slot-positions of the top-left element of b on a.

# Get range of positions (rows and cols) where we have zeros
target_indices = np.argwhere(a==0)
minmax = np.array([target_indices.min(axis=0), target_indices.max(axis=0)])
# Define max position (index) of topleft element of b on a
maxpos = np.dot(np.array([-1,1]), minmax) + minmax[0,:] - (np.array(b.shape) -1)
# Define min position (index) of topleft element of b on a
minpos = minmax[0,:]

Make a list of Top Left Corner Positions of b on a

Th function get_rand_topleftpos() takes in minpos and maxpos for rows and columns on a that define possible slot-positions, and returns a randomly selected valid slot-position for size=1. I have used a size=20 to create quite a few valid random slot-positions and then select only the unique positions so we could then see them as images. If you need just one slot-position at a time, select size=1.

def get_rand_topleftpos(minpos, maxpos, size=1):
    rowpos = np.random.randint(minpos[0], high=maxpos[0] + 1, size=size)
    colpos = np.random.randint(minpos[1], high=maxpos[1] + 1, size=size)
    pos = np.vstack([rowpos, colpos]).T
    return (rowpos, colpos, pos)

# Make a few valid positions where the array b could be placed
rowpos, colpos, pos = get_rand_topleftpos(minpos, maxpos, size=20)
# Select the Unique combinations so we could visualize them only
pos = np.unique(pos, axis=0)

Place b on a and Make Figures

We make a custom defined function fill_a_with_b() to fill a with b at a certain postion on a. This position will accept the top-left cell of b.

def fill_a_with_b(a, b, pos = [0,0]):
    aa = a.copy()
    aa[slice(pos[0], pos[0] + b.shape[0]), 
       slice(pos[1], pos[1] + b.shape[1])] = b.copy()
    return aa

# Make a few figures with randomly picked position 
# for topleft position of b on a

if pos.shape[0]>6:
    nrows, ncols = int(np.ceil(pos.shape[0]/6)), 6
else:
    nrows, ncols = 1, pos.shape[0]
fig, axs = plt.subplots(nrows = nrows, 
                        ncols = ncols, 
                        figsize=(2.5*ncols,2.5*nrows))
for i, ax in enumerate(axs.flatten()):    
    if i<pos.shape[0]:
        aa = fill_a_with_b(a, b, pos[i,:]) 

        sns.heatmap(aa, 
                    vmin=np.min(aa), 
                    vmax=np.max(aa), 
                    annot=True, 
                    cbar=False,
                    square=True,
                    cmap = 'YlGnBu_r', 
                    ax = ax
                   ); 

        ax.set_title('TopLeftPos: {}'.format(tuple(pos[i,:])), 
                     fontsize=9);
    else:
        ax.axis('off')

plt.tight_layout()        
plt.show()

Results

With array a defined as:

shape = (5,5)
a = np.zeros(shape)
a[:,-1] = np.arange(shape[0]) + 10
a[-1,:] = np.arange(shape[1]) + 10

With array a defined as:

shape = (6,5)
a = np.zeros(shape)
a[:,0] = np.arange(shape[0]) + 10
a[:,-1] = np.arange(shape[0]) + 10
a[-1,:] = np.arange(shape[1]) + 10