I realized that I hadn't correctly formatted the predictions made by the network with dlarray. Below is the working regression layer based on Huber loss in case anyone needs it.
classdef huberRegressionLayer < nnet.layer.RegressionLayer ...
& nnet.layer.Acceleratable % (Optional)
methods
function layer = huberRegressionLayer(name)
% layer = huberRegressionLayer(name) creates a
% robust regression layer based on huber loss
% and specifies the layer name.
% Set layer name.
layer.Name = name;
% Set layer description.
layer.Description = 'Huber loss';
end
function loss = forwardLoss(layer,Y,T)
% loss = forwardLoss(layer, Y, T) returns the Huber loss between
% the predictions Y and the training targets T.
dlY = dlarray(Y,'CB');
loss = huber(dlY,T,"TransitionPoint",1);
end
end
end