Is linespace supposed to be single precision only?

[cianciosa]

I'm trying to create uniform spaced array of angles from 0 to Pi but I'm running into a precision issues. The arrays that I'm creating are being reported as float32 and float64 respectively.

mlx::core::set_default_device(mlx::core::Device::cpu);

mlx::core::array float_angles  = mlx::core::linspace(0, std::numbers::pi_v<float>,  5, mlx::core::float32);
mlx::core::array double_angles = mlx::core::linspace(0, std::numbers::pi_v<double>, 5, mlx::core::float64);
float_angles.eval();
double_angles.eval();

std::cout << float_angles << std::endl << double_angles << std::end;

This produces the output

array([0, 0.785398, 1.5708, 2.35619, 3.14159], dtype=float32)
array([0, 0.785398, 1.5708, 2.35619, 3.14159], dtype=float64)

However, when I compare these values against what I calculate manually, it appears the array that reports its type as float64 actually contains values which are float32.

std::array<float, 5> cpp_float_angles;
std::array<double, 5> cpp_double_angles;
cpp_float_angles[0] = 0;
cpp_double_angles[0] = 0;

const float dtheta_float = std::numbers::pi_v<float>/4;
const double dtheta_double = std::numbers::pi_v<double>/4;
for (size_t i = 1; i < 5; i++) {
    cpp_float_angles[i] = i*dtheta_float;
    cpp_double_angles[i] = i*dtheta_double;
}
    
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10);
for (size_t i = 0; i < 5; i++) {
    std::cout << *(float_angles.data<float> () + i) - cpp_float_angles[i] << " "
              << *(double_angles.data<double> () + i) - cpp_double_angles[i] << " "
              << *(double_angles.data<double> () + i) - cpp_float_angles[i] << " "
              << cpp_float_angles[i] - cpp_double_angles[i] << std::endl;
}

I get the result showing the float64 arrays has values identical to the single precision values.

0 0 0 0
0 2.1855695031547384e-08 0 2.1855695031547384e-08
0 4.3711390063094768e-08 0 4.3711390063094768e-08
0 5.9624403192515274e-09 0 5.9624403192515274e-09
0 8.7422780126189537e-08 0 8.7422780126189537e-08

[awni]

That's a bug, thanks for flagging. I will send a patch in a moment.

[cianciosa]

The patch fixes the issue with linspace so I started checking some operations on that array now. I wouldn't expect these to match exactly in all cases since the implementations may be different between the library MLX is using the libc++. But there appears to be times when results revert back to single precision.

#define OP abs

mlx::core::set_default_device(mlx::core::Device::cpu);

mlx::core::array float_angles  = mlx::core::linspace(0, std::numbers::pi_v<float>,  5, mlx::core::float32);
mlx::core::array double_angles = mlx::core::linspace(0, std::numbers::pi_v<double>, 5, mlx::core::float64);
float_angles.eval();
double_angles.eval();

mlx::core::array op_float = mlx::core::OP(float_angles);
mlx::core::array op_double = mlx::core::OP(double_angles);
op_float.eval();
op_double.eval();

std::cout << op_float << std::endl << op_double << std::endl << std::endl;

std::array<float, 5> cpp_float_ops;
std::array<double, 5> cpp_double_ops;
cpp_float_ops[0] = std::OP<float> (0);
cpp_double_ops[0] = std::OP<double> (0);

const float dtheta_float = std::numbers::pi_v<float>/4;
const double dtheta_double = std::numbers::pi_v<double>/4;
for (size_t i = 1; i < 5; i++) {
    cpp_float_ops[i] = std::OP<float> (i*dtheta_float);
    cpp_double_ops[i] = std::OP<double> (i*dtheta_double);
}
    
std::cout << std::setprecision(std::numeric_limits<double>::max_digits10);
for (size_t i = 0; i < 5; i++) {
    std::cout << *(op_float.data<float> () + i) - *(op_double.data<double> () + i) << " "
              << *(op_float.data<float> () + i) - cpp_float_ops[i] << " "
              << *(op_double.data<double> () + i) - cpp_float_ops[i] << " "
              << cpp_float_ops[i] - cpp_double_ops[i] << std::endl;
}

Here are the results for the operations that show some anomalies

exp

array([1, 2.19328, 4.81048, 10.5507, 23.1407], dtype=float32)
array([1, 2.19328, 4.81048, 10.5507, 23.1407], dtype=float64)

0 0 0 0 0
0 0 -6.9124856061364426e-08 0 -6.9124856061364426e-08
0 -4.76837158203125e-07 -1.2419044903566601e-07 -4.76837158203125e-07 3.5264670916745899e-07
0 0 -4.4656745856741509e-08 0 -4.4656745856741509e-08
0 -3.814697265625e-06 -8.7557711836439012e-07 -3.814697265625e-06 2.9391201472606099e-06

Note for this case, I can verify this is single precision in an alternative method. If I take the exponential of a large number there is a point where single precision should overflow but double would not.

 const double d_arg = 600;
 const float f_arg = 600;
 const double d_result = std::exp<double> (d_arg);
 const float f_result = std::exp<float> (f_arg);
 std::cout << d_result << " " << f_result << std::endl;
 std::cout << std::numeric_limits<float>::max() << std::endl;
 std::cout << std::numeric_limits<double>::max() << std::endl;

 mlx::core::array test_float = mlx::core::exp(mlx::core::array(600, mlx::core::float64));
 mlx::core::array test_double = mlx::core::exp(mlx::core::array(600, mlx::core::float32));
 std::cout << test_float << test_double << std::endl;

It is expected that single precision would be inf but double precision should not.

3.7730203009299397e+260 inf
3.4028234663852886e+38
1.7976931348623157e+308
array(inf, dtype=float64)array(inf, dtype=float32)

sin

array([0, 0.707107, 1, 0.707107, -8.74228e-08], dtype=float32)
array([0, 0.707107, 1, 0.707107, -8.74228e-08], dtype=float64)

0 0 0 0 0
0 5.9604644775390625e-08 4.7503027733597492e-08 5.9604644775390625e-08 -1.2101617041793133e-08
0 0 0 0 0
0 0 -1.2101617152815436e-08 0 -1.2101617152815436e-08
0 0 -8.7422776695940543e-08 0 -8.7422776695940543e-08

cos

array([1, 0.707107, -4.37114e-08, -0.707107, -1], dtype=float32)
array([1, 0.707107, -4.37114e-08, -0.707107, -1], dtype=float64)

0 0 0 0 0
0 0 -1.2101617152815436e-08 0 -1.2101617152815436e-08
0 0 -4.3711388347970271e-08 0 -4.3711388347970271e-08
0 0 1.2101617041793133e-08 0 1.2101617041793133e-08
0 0 0 0 0

erf

array([0, 0.733311, 0.973679, 0.999138, 0.999991], dtype=float32)
array([-0, 0.733311, 0.973679, 0.999138, 0.999991], dtype=float64)

0 0 -0 -0 0
0 -5.9604644775390625e-08 -7.0451928957382393e-08 -5.9604644775390625e-08 -1.0847284181991768e-08
0 5.9604644775390625e-08 8.1021442444040304e-08 5.9604644775390625e-08 2.1416797668649679e-08
0 -5.9604644775390625e-08 -8.3415520490959238e-08 -5.9604644775390625e-08 -2.3810875715568613e-08
0 0 -4.9457038464595371e-09 0 -4.9457038464595371e-09