Suggestion: Allow for some tolerance in numerical prescion in computation of x^T H x test calculations

[robbietuk]

A suggestion based upon a failing test thread in stir-users email list.

This objective function concavity test was found to fail:

Error : 0 is larger than 0,
FAIL: PoissonLLListModeData: Computation of x^T H x = 0.000000 > 0 (Hessian) and is therefore NOT concave
 >target image max=0.999849
 >target image min=1.16359e-05
 >output image max=0
 >output image min=0
Writing diagnostic files to PoissonLLListModeData_concavity_out.hv, *target.hv

that comes from here:

STIR/src/include/stir/recon_buildblock/test/ObjectiveFunctionTests.h

Lines 228 to 244 in a008e73

	// test for a CONCAVE function (after multiplying with mult_factor of course)
	if (this->check_if_less(my_sum, 0))
	{
	return Succeeded::yes;
	}
	else
	{
	// print to console the FAILED configuration
	std::cerr << "FAIL: " + test_name + ": Computation of x^T H x = " + std::to_string(my_sum)
	<< " > 0 (Hessian) and is therefore NOT concave"
	<< "\n >target image max=" << target.find_max() << "\n >target image min=" << target.find_min()
	<< "\n >output image max=" << output->find_max() << "\n >output image min=" << output->find_min() << '\n';
	std::cerr << "Writing diagnostic files to " << test_name + "_concavity_out.hv, *target.hv\n";
	write_to_file(test_name + "_concavity_out.hv", *output);
	write_to_file(test_name + "_target.hv", target);
	return Succeeded::no;
	}

This isn't found in CI but the test might be too strict. I suggest the x^T H x calculation is compared to a small negative number, similar to that in test_priors

STIR/src/recon_test/test_priors.cxx

Lines 247 to 268 in a008e73

	/// Compute x \cdot (H x)
	const double my_sum = std::inner_product(input->begin_all(), input->end_all(), output->begin_all(), double(0));
	/// Compute x \cdot x
	const double my_norm2 = std::inner_product(input->begin_all(), input->end_all(), input->begin_all(), double(0));
	// test for a CONVEX function: 0 < my_sum, but we allow for some numerical error
	if (this->check_if_less(-my_norm2 * 2E-4, my_sum))
	{
	return true;
	}
	else
	{
	// print to console the FAILED configuration
	info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " < 0 and is therefore NOT convex" + "\ntest_name="
	+ test_name + "\nbeta=" + std::to_string(beta) + "\ninput_multiplication=" + std::to_string(input_multiplication)
	+ "\ninput_addition=" + std::to_string(input_addition) + "\ncurrent_image_multiplication="
	+ std::to_string(current_image_multiplication) + "\ncurrent_image_addition=" + std::to_string(current_image_addition)
	+ "\n >input image max=" + std::to_string(input->find_max()) + "\n >input image min="
	+ std::to_string(input->find_min()) + "\n >input image norm^2=" + std::to_string(my_norm2) + "\n >target image max="
	+ std::to_string(target.find_max()) + "\n >target image min=" + std::to_string(target.find_min()));
	return false;
	}

[KrisThielemans]

This failure occurs on MacOS CI, which is therefore disabled... See #1491 for more info and some other cases where it fails.

It does appear that the test should be relaxed

STIR/src/include/stir/recon_buildblock/test/ObjectiveFunctionTests.h

Line 229 in a008e73

if (this->check_if_less(my_sum, 0))

tests for "strict" concavity, which is not guaranteed. I wonder if using check_if_less_or_equal (which actually we don't have in RunTests) would resolve this.

On the other hand, I'm surprised that this tests gives 0 for x^T H x. I had a look at the construction of target, which seems to be here:

STIR/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx

Lines 147 to 154 in a008e73

	write_to_file("target.hv", *density_sptr);
	// fill with random numbers between 0 and 1
	typedef boost::mt19937 base_generator_type;
	// initialize by reproducible seed
	static base_generator_type generator(boost::uint32_t(42));
	static boost::uniform_01<base_generator_type> random01(generator);
	for (target_type::full_iterator iter = density_sptr->begin_all(); iter != density_sptr->end_all(); ++iter)
	*iter = static_cast<float>(random01());

I note that the writing is in the wrong place (should be after the filling), so we currently cannot check what's in there, but according to the min/max in the diagnostics, target is (strictly) positive. That implies that H x has to be non-negative, and hence x^T H x > 0.

In test_priors, we needed to allow a small negative number because H contains negative numbers. However, for the Poisson LL, this is not the case.

So, I don't understand what's going on, and am therefore reluctant to just let the test pass...

[robbietuk]

That makes sense. I'll close this because #1491 probably covers it.