What is the fastest method to searsch for a part of a matrix inside of a large matrix?

Assume that you have this matrix

X = randn(100, 100);

And you cut out this part

% Windows 95
m = randi(95, 1)
n = randi(95, 1)
x = X(m:m+5, n:n+5);

Question:

In the real world, it's very naive to think that it's possible to have a slinding window for scanning where x can fit in X and check the smallest error.

smallest = inf;
for i = 1:95
  for j = 1:95
    % Extract candidate 
    candidate = X(i:i+5, j:j+5);
    
    % Compare 
    difference = abs(candidate - x);
    
    % Find error 
    error = sum(difference(:));
    if(error < smallest)
      smallest = error;
      index_i = i;
      index_j = j;
    end 
  end
end

fprintf('The smallest error was %f and it was located at %i, %i', smallest, index_i, index_j);

A small run of this code would be:

>> m = randi(95, 1)
  m = 21
>> n = randi(95, 1)
  n = 71
>> x = X(m:m+5, n:n+5);
>> The smallest error was 0.000000 and it was located at 21, 71

I want to find candidates of x that might fit into X, but I don't know which algorithm I should use.

The purpose of this question is that I'm going to do very basic landmark detection for images. Like this image below.

Normally, Neural Networks are used for landmark detection, but I don't have so much data so I could train a neural network. So I need to rely on other methods. Which method should I use?

This would be done with a convolution, as I think you are aware because you used the convolution tag. In Matlab, that's best done using normxcorr2. Here's an example of some code:

% Set dimensions
Nx = 1000;
Ny = 1000;
SmallNx = 20;

% Generate matrix
X = rand([Nx Ny]);

% Region of interest
xi = randi(Nx - SmallNx);
yi = randi(Ny - SmallNx);
SmallX = X(xi:xi+SmallNx, yi:yi+SmallNx);

% Convolution
ConvX = normxcorr2(SmallX, X);

% Find convolution peak, correct for padding
[xi_est, yi_est] = find(ConvX == max(max(ConvX)));
xi_est = xi_est - SmallNx;
yi_est = yi_est - SmallNx;

For a more intuitive example, you could use an image:

% Load image
load Mandrill;
Nx = size(X, 1);
Ny = size(X, 2);
SmallNx = 100;

% Region of interest
xi = randi(Nx - SmallNx);
yi = randi(Ny - SmallNx);
SmallX = X(xi:xi+SmallNx, yi:yi+SmallNx);

% Convolution
ConvX = normxcorr2(SmallX, X);

% Find convolution peak, correct for padding
[xi_est, yi_est] = find(ConvX == max(max(ConvX)));
xi_est = xi_est - SmallNx;
yi_est = yi_est - SmallNx;

figure; 
subplot(1,2,1);
imagesc(X);
colormap gray;
hold on
rectangle('Position', [yi_est xi_est SmallNx SmallNx], 'EdgeColor','r')
subplot(1,2,2);
imagesc(SmallX);
colormap gray;