Different behavior between vector-wise normalize and base vector normalize

Summary

The implementation of the normalize-related functions for MatrixBase preserves the original values when the squared norm is zero, see: https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/Core/Dot.h?ref_type=heads#L92

However, the same functions produce NaNs when applied on top of colwise/rowwise methods, see https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/Core/VectorwiseOp.h?ref_type=heads#L633.

This discrepancy is not explicitly documented and may lead to bugs if users expect the vector-wise operations to always match MatrixBase behavior.

I am not very familiar with how the expression templates in Eigen work, but I would expect VectorwiseOp::normalized() to return an expression that refers to the block normalize methods. I assume there may be good reasons to implement it the current way, using Replicate and cwise division by the norm, but I would be happy to know from Eigen experts if there is a way to alleviate this issue in the future or at least to add some warning in the documentation regarding the matters.

Environment

Should be reproducible regardless of the environment. We used version 3.4.0 of Eigen when we encountered issues, but I believe the issue persists in actual development code.

Minimal Example

// These two functions produce different results:
void normalize1(Eigen::MatrixXf& matrix)
{
    for (auto&& col : matrix.colwise()) col.normalize();
}
void normalize2(Eigen::MatrixXf& matrix)
{
    matrix.colwise().normalize();
}

Steps to reproduce

Create a matrix that has columns of zero norm, run the code containing the two functions above, compare the results.

What is the current bug behavior?

The normalize behavior differs between the for loop and the vector-wise method.

What is the expected correct behavior?

The normalize behavior is the same (or almost the same, given floating point errors) as the vector-wise method.

Anything else that might help

I envision that an ultimate fix could be one of the following:

• Implement a vector-wise expression type that is a view over normalized columns/rows from the block normalized() method (not sure it is feasible in current Eigen philosophy)
• Add an unary expression that replaces zeros with ones for the divisor (performance considerations)

I understand, however, that the behavior change versus the old version may be not appreciated, but in my opinion, the existing discrepancy may be more severe. Highlighting the differences in the docs could be helpful as well.