Is MSE loss a valid ELBO loss to measure? I am learning from an example given by TensorFlow document, https://www.tensorflow.org/tutorials/generative/cvae#define_the_loss_function_and_the_optimizer:

VAEs train by maximizing the evidence lower bound (ELBO) on the
marginal log-likelihood.
In practice, optimize the single sample Monte Carlo estimate of this
expectation: logp(x|z) + logp(z) - logq(z|x).

The loss function was implemented as:
def log_normal_pdf(sample, mean, logvar, raxis=1):
  log2pi = tf.math.log(2. * np.pi)
  return tf.reduce_sum(
      -.5 * ((sample - mean) ** 2. * tf.exp(-logvar) + logvar + log2pi),
      axis=raxis)


def compute_loss(model, x):
  mean, logvar = model.encode(x)
  z = model.reparameterize(mean, logvar)
  x_logit = model.decode(z)
  cross_ent = tf.nn.sigmoid_cross_entropy_with_logits(logits=x_logit, labels=x)
  logpx_z = -tf.reduce_sum(cross_ent, axis=[1, 2, 3])
  logpz = log_normal_pdf(z, 0., 0.)
  logqz_x = log_normal_pdf(z, mean, logvar)
  return -tf.reduce_mean(logpx_z + logpz - logqz_x)

Since this example used MINIST dataset, x can be normalized to [0, 1] and sigmoid_cross_entropy_with_logits was used here.
My question:
Another example used MSE loss (as follow), is MSE loss a valid ELBO loss to measure p(x|z)? Can we use other loss functions as a reconstruction loss in VAE, such as Huber loss (https://en.wikipedia.org/wiki/Huber_loss)?
https://www.tensorflow.org/guide/keras/custom_layers_and_models#putting_it_all_together_an_end-to-end_example
    # Iterate over the batches of the dataset.
    for step, x_batch_train in enumerate(train_dataset):
        with tf.GradientTape() as tape:
            reconstructed = vae(x_batch_train)
            # Compute reconstruction loss
            loss = mse_loss_fn(x_batch_train, reconstructed)
            loss += sum(vae.losses)  # Add KLD regularization loss


As you point out, the ELBO is approximated by $\log p(x|z) + \log p(z) - \log q(z|x)$ with reparameterized samples from $q(z|x)$.
The first term, $\log p(x|z)$, is, in expectation over $q(z|x)$, often called the likelihood loss (reconstruction loss is a synonym in the context of VAEs). Such a loss, but without the expectation over $q(z|x)$, appears also in neural networks without latent variables $z$, in the form of $\log p(y|x)$, e.g. in classification/regression.
We have freedom over the choice of the modeled distribution $p(x|z)$, and depending on that distribution get a different loss function.
The likelihood loss for a Bernoulli output distribution (per pixel) is sigmoid_cross_entropy_with_logits and for a Gaussian output distribution (per pixel) is the MSE loss.
The latter is because (up to an additive and a multiplicative constant) the quadratic error equals the log of the Gaussian likelihood/density.
Similarly for other loss functions, also for the Huber loss, you can find a distribution $p(x|z)$ that corresponds to it. It is therefore completely reasonable to use any such loss functions.
However, the loss function should fit the output domain. If it's discrete, you shouldn't use a continuous loss function and vice versa. But with the Huber loss for continuous values, you are on the right track.
It's also not very elegant to use an infinite range output distribution / corresponding loss if the output is on a finite range like pixel values, but since all models are wrong to some extend anyway, no need to bother I guess, if your library doesn't offer something for finite ranges. Note that if your values are between 0 and 255, it's always possible to normalize them if that is convenient.
