function betaDivPerElem = betaDivergencePerElement(x, y, beta)

switch beta
    case 0
        betaDivPerElem = x./y - log(x./y) - 1;

    case 1
        betaDivPerElem = x .* log(x./y) + x - y;

    otherwise
        betaDivPerElem = (x.^beta) / (beta*(beta-1)) ...
                         + (y.^beta) / beta ...
                         - (x .* (y.^(beta-1))) / (beta-1);
end