def METHOD_NAME(A: np.ndarray, B: np.ndarray, **kwargs: Any) -> np.ndarray:
    """Non-negative least squares.
    Given two matrices A and B, find a non-negative matrix X
    that minimizes the sum squared error::
        err(X) = sum_i,j ((AX)[i,j] - B[i, j])^2
    Parameters
    ----------
    A : np.ndarray [shape=(m, n)]
        The basis matrix
    B : np.ndarray [shape=(..., m, N)]
        The target array.  Additional leading dimensions are supported.
    **kwargs
        Additional keyword arguments to `scipy.optimize.fmin_l_bfgs_b`
    Returns
    -------
    X : np.ndarray [shape=(..., n, N), non-negative]
        A minimizing solution to ``|AX - B|^2``
    See Also
    --------
    scipy.optimize.nnls
    scipy.optimize.fmin_l_bfgs_b
    Examples
    --------
    Approximate a magnitude spectrum from its mel spectrogram
    >>> y, sr = librosa.load(librosa.ex('trumpet'), duration=3)
    >>> S = np.abs(librosa.stft(y, n_fft=2048))
    >>> M = librosa.feature.melspectrogram(S=S, sr=sr, power=1)
    >>> mel_basis = librosa.filters.mel(sr=sr, n_fft=2048, n_mels=M.shape[0])
    >>> S_recover = librosa.util.nnls(mel_basis, M)
    Plot the results
    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(nrows=3, sharex=True, sharey=True)
    >>> librosa.display.specshow(librosa.amplitude_to_db(S, ref=np.max),
    ...                          y_axis='log', x_axis='time', ax=ax[2])
    >>> ax[2].set(title='Original spectrogram (1025 bins)')
    >>> ax[2].label_outer()
    >>> librosa.display.specshow(librosa.amplitude_to_db(M, ref=np.max),
    ...                          y_axis='mel', x_axis='time', ax=ax[0])
    >>> ax[0].set(title='Mel spectrogram (128 bins)')
    >>> ax[0].label_outer()
    >>> img = librosa.display.specshow(librosa.amplitude_to_db(S_recover, ref=np.max(S)),
    ...                          y_axis='log', x_axis='time', ax=ax[1])
    >>> ax[1].set(title='Reconstructed spectrogram (1025 bins)')
    >>> ax[1].label_outer()
    >>> fig.colorbar(img, ax=ax, format="%+2.0f dB")
    """
    # If B is a single vector, punt up to the scipy method
    if B.ndim == 1:
        return scipy.optimize.METHOD_NAME(A, B)[0]  # type: ignore
    n_columns = int(MAX_MEM_BLOCK // (np.prod(B.shape[:-1]) * A.itemsize))
    n_columns = max(n_columns, 1)
    # Process in blocks:
    if B.shape[-1] <= n_columns:
        return _nnls_lbfgs_block(A, B, **kwargs).astype(A.dtype)
    x: np.ndarray
    x = np.einsum("fm,...mt->...ft", np.linalg.pinv(A), B, optimize=True)
    np.clip(x, 0, None, out=x)
    x_init = x
    for bl_s in range(0, x.shape[-1], n_columns):
        bl_t = min(bl_s + n_columns, B.shape[-1])
        x[..., bl_s:bl_t] = _nnls_lbfgs_block(
            A, B[..., bl_s:bl_t], x_init=x_init[..., bl_s:bl_t], **kwargs
        )
    return x