spectrum: computing lambda fails when reducedPrefactor=inf

[vasdommes]

When reducedPrefactor has a pole at x=0 and spectrum finds zero at x=0, then spectrum crashes when writing lambda to file:

GNU MP: Cannot allocate memory (size=16140901064495857912)
[exp-9-55:3563804] *** Process received signal ***
[exp-9-55:3563804] Signal: Aborted (6)
[exp-9-55:3563804] Signal code:  (-6)
[exp-9-55:3563804] [ 0] /lib64/libpthread.so.0(+0x12cf0)[0x15554ae9acf0]
[exp-9-55:3563804] [ 1] /lib64/libc.so.6(gsignal+0x10f)[0x15554ab11acf]
[exp-9-55:3563804] [ 2] /lib64/libc.so.6(abort+0x127)[0x15554aae4ea5]
[exp-9-55:3563804] [ 3] /home/vdommes/install/lib/libgmp.so.10(+0xdb22)[0x155554eb6b22]
[exp-9-55:3563804] [ 4] /home/vdommes/install/lib/libgmp.so.10(__gmp_tmp_reentrant_alloc+0x21)[0x155554eb76c1]
[exp-9-55:3563804] [ 5] /home/vdommes/install/lib/libgmp.so.10(__gmpf_get_str+0x791)[0x155554eb8f21]
[exp-9-55:3563804] [ 6] /home/vdommes/install/lib/libgmp.so.10(__gmp_doprnt_mpf2+0x67)[0x155554f03617]
[exp-9-55:3563804] [ 7] /home/vdommes/install/lib/libgmpxx.so.4(_ZlsRSoPK12__mpf_struct+0xb1)[0x155555122a21]
[exp-9-55:3563804] [ 8] /home/vdommes/install/lib64/libEl.so.0(_ZlsIA1_12__mpf_structS1_ERSoS2_RK10__gmp_exprIT_T0_E+0x36)[0x1555521c62d6]
[exp-9-55:3563804] [ 9] /home/vdommes/install/lib64/libEl.so.0(_ZN2EllsERSoRKNS_8BigFloatE+0x23)[0x1555521c46d0]
[exp-9-55:3563804] [10] /home/vdommes/install/sdpb/master/bin/spectrum[0x47c476]
[exp-9-55:3563804] [11] /home/vdommes/install/sdpb/master/bin/spectrum[0x47b2bc]
[exp-9-55:3563804] [12] /home/vdommes/install/sdpb/master/bin/spectrum[0x420d7c]
[exp-9-55:3563804] [13] /lib64/libc.so.6(__libc_start_main+0xe5)[0x15554aafdd85]
[exp-9-55:3563804] [14] /home/vdommes/install/sdpb/master/bin/spectrum[0x424aae]

Example input:
pmp_info.json
out.zip
Command line:
spectrum --lambda 1 --precision 768 --threshold 1e-10 --pmpInfo pmp_info.json --solution out --output spectrum.json

[vasdommes]

Explanation

The error originates from this line:

sdpb/src/spectrum/compute_spectrum/compute_lambda.hxx

Lines 207 to 213 in f02430b

	// Set lambda = 1/sqrt(\chi) * v_{j,\tau}
	// With this definition, lambda does not change
	// if one adds reducedPrefactor != prefactor to PMP.json
	// NB: this is different from Python script and from (A.8) definition!
	lambda *= 1
	/ El::Sqrt(to_BigFloat(pvm_info.reduced_prefactor.evaluate(
	to_Boost_Float(zeros.back().zero))));

Here chi = reduced_prefactor.evaluate(0) returns inf due to a pole at x=0.

In theory, dividing by sqrt(inf) yields zero. But GMP cannot handle infinities, so to_BigFloat(inf) is an invalid GMP number (with a huge invalid exponent). The error propagates silently, and we get lambda = inf instead of lambda = 0.
GMP is trying to allocate an absurd amount of memory to print a huge invalid number, and crashes.
Note that if wrap lambda with to_Boost_Float() here:

sdpb/src/spectrum/write_spectrum/write_file.cxx

Line 55 in 4a38186

<< zeros_iterator->zeros.at(zero_index).lambda(row, 0)

it will print inf instead of crashing.

[vasdommes]

Note also that this zero should be excluded because x = 0 is not a minimum. But this line:

sdpb/src/spectrum/compute_spectrum/find_zeros.cxx

Line 259 in fec8e93

is_zero = y / y_right < threshold;

works incorrectly and sets is_zero = true when y = inf due to silent GMP errors (in the test case above, for block_index=12 it essentially computed to_BigFloat(inf) / 0.00588129 = 9.21738e-18).

y is computed as a det(interpolation_poly_matrix(x) * reduced_prefactor(x)).

In principle, it is also possible to have infinite prefactor and zero determinant, so that their product is NaN.

[vasdommes]

Possible fixes

MPFR -> GMP

First of all, in to_Boost_Float() we should always check that value is a valid number, e.g. via ASSERT(isfinite(value)).

spectrum

Then, we have to decide what to do in spectrum when chi = inf. NB: this happens not only in compute_lambda.hxx, but also in find_zeros.cxx, when we evaluate reduced_prefactor at all zeros found by MPSolve and also at x = 0.

find_zeros (check minima)

sdpb/src/spectrum/compute_spectrum/find_zeros.cxx

Lines 246 to 297 in fec8e93

	Scoped_Timer check_minima_timer(timers, "check_minima");
	for(size_t i = 0; i < minima.size(); ++i)
	{
	const auto &x = minima.at(i);
	const auto y = eval(x);
	bool is_zero = false;
	if(i == 0)
	{
	if(minima.size() > 1)
	{
	const auto x_right = get_midpoint(x, minima.at(i + 1));
	const auto y_right = eval(x_right);
	is_zero = y / y_right < threshold;
	}
	else
	{
	// This is a case of single minimum.
	// TODO: which points should we choose for comparison?
	// It's not obvious, choosing x/2 is not justified well.
	auto x_other = x / 2;
	if(x_other == El::BigFloat(0))
	{
	// Special case: if x=0, we take the first nonzero sample point.
	x_other = pvm.sample_points.at(0);
	if(x_other == El::BigFloat(0))
	x_other = pvm.sample_points.at(1);
	}
	ASSERT(x_other > 0);
	const auto y_other = eval(x_other);
	is_zero = y / y_other < threshold;
	}
	}
	else if(i + 1 == minima.size())
	{
	const auto x_left = get_midpoint(x, minima.at(i - 1));
	const auto y_left = eval(x_left);
	is_zero = y / y_left < threshold;
	}
	else
	{
	const auto x_left = get_midpoint(x, minima.at(i - 1));
	const auto y_left = eval(x_left);
	const auto x_right = get_midpoint(x, minima.at(i + 1));
	const auto y_right = eval(x_right);
	is_zero = y * y / y_left / y_right < threshold * threshold;
	}
	if(is_zero)
	zeros.emplace_back(x);
	}

In find_zeros.cxx, we compute det(c - By) at all roots found by MPSolve and also at x = 0.

By comparing values of det(c - By) at neighboring points, we decide which of them are minima that are deep enough (below --threshold).

The question is, what should we do if some values are inf or NaN?
For example, what if we have a sequence [inf, 1e-10, NaN, 100]? Is 1e-10 a valid minimum? Shall we throw an error at NaN, or drop it silently?
(In practice, we expect only inf at x = 0 and ordinary numbers at all other locations)

compute_lambda

Currently, we compute contribution error_matrix first and then multiply lambda by 1/sqrt(\chi):

sdpb/src/spectrum/compute_spectrum/compute_lambda.hxx

Lines 200 to 213 in fec8e93

	error_matrix(index, row_column)
	-= interpolation(index, zero_index) * lambda(row)
	* lambda(column) * (row == column ? 1 : 2);
	}
	++row_column;
	}
	// Set lambda = 1/sqrt(\chi) * v_{j,\tau}
	// With this definition, lambda does not change
	// if one adds reducedPrefactor != prefactor to PMP.json
	// NB: this is different from Python script and from (A.8) definition!
	lambda *= 1
	/ El::Sqrt(to_BigFloat(pvm_info.reduced_prefactor.evaluate(
	to_Boost_Float(zeros.back().zero))));

If lambda becomes zero, shall we exit early, without updating error_matrix?
Note that in the example above max_eigenvalue is small (~ dualityGap = 1e-50), so it does not really matter.
If it is big, then not updating error_matrix will make error large, indicating possible issues.
Maybe it is worth printing a warning if error > threshold.