This is closely related to Blind Deconvolution.
The only difference is we limit our self to a very specific type of blur kernels.
The nice thing about the Gaussian Kernel is being defined by single parameter - The Standard Deviation of the kernel.
The less nice thing is the connection isn't linear.
Optimization Problem
Let's define a classic non linear model for this problem:
$$ \arg \min_{\sigma} \frac{1}{2} {\left\| A \left( \sigma \right) x - b \right\|}_{2}^{2} $$
Where $ A \left( \sigma \right) $ is the convolution matrix generated by a Gaussian Kernel parameterized by $ \sigma $, $ x $ is the original image ($ A $ in your question) in a vector shape, and $ b $ is the blurred image ($ B $ in your question).
This is a classical Non Linear Least Squares problem which can be solved by MATLAB using lsqnonlin()
.
Code Sample
This is the main part of the code:
mA = rand([numRows, numCols]);
gaussianKernelStd = kernelStdLowerBound + ((kernelStdUpperBound - kernelStdLowerBound) * rand(1));
mB = ApplyGaussianBlur(mA, gaussianKernelStd, STD_TO_RADIUS_FACTOR);
% Objective Functions
hObjFun = @(kernelStd) reshape(ApplyGaussianBlur(mA, kernelStd, STD_TO_RADIUS_FACTOR) - mB, [numPx, 1]);
estKernelEst = lsqnonlin(hObjFun, initKernelStd, kernelStdLowerBound, kernelStdUpperBound, sSolverOptions);
Results

As can be seen, the estimation is almost perfect.
On large images it might take time (Using some tricks of the Gaussian Filter in the Fourier Domain one could do that there with major speed up), but still, it is not free.
The full code is available on my StackExchange Signal Processing Q49121 GitHub Repository (Look at the SignalProcessing\Q49121
folder).