Блог

Я написал класс-оболочку C для некоторых функций в LAPACK. Чтобы протестировать класс, я использую расширение Python C, куда я вызываю numpy , и выполняю те же операции, и сравниваю результаты, принимая разницу

Например, для обратной матрицы я генерирую случайную матрицу на C , затем передаю ее в виде строки (со многими, многими цифрами, например, 30 цифр) на терминал Python с использованием PyRun_SimpleString и присваиваю матрице значение numpy.matrix(...,dtype=numpy.double) (или numpy.complex128 ). Затем я использую numpy.linalg.inv() для вычисления обратной той же матрицы. Наконец, я беру разницу между numpy результатом и моим результатом и использую numpy.isclose с определенным относительным допуском, чтобы увидеть, достаточно ли близки результаты.

Проблема: проблема в том, что когда я использую C с плавающей точкой, относительная точность, которую мне нужно сравнить, составляет около 1e-2 !!! И все же с этой относительной точностью я получаю некоторые статистические сбои (с низкой вероятностью).

Удвоения в порядке… Я могу это сделать 1e-10 , и это статистически безопасно.

Хотя я знаю, что значения с плавающей запятой имеют внутреннюю битовую точность около 1e-6 , мне интересно, почему я должен опускаться так низко 1e-2 , чтобы иметь возможность сравнивать результаты, и это все еще терпит неудачу несколько раз!

Итак, опускаясь так низко 1e-2 , я задался вопросом, не думаю ли я обо всем этом неправильно. Что-то не так с моим подходом?

Пожалуйста, запросите более подробную информацию, если вам это нужно.

Обновление 1: Эрик запросил пример вызовов Python. Вот пример:

 //create my matrices
Matrix<T> mat_d = RandomMatrix<T>(...);
auto mat_d_i = mat_d.getInverse();

//I store everything in the dict 'data'
PyRun_SimpleString(std::string("data={}").c_str());

//original matrix
//mat_d.asString(...) will return in the format [[1,2],[3,4]], where 32 is 32 digits per number
PyRun_SimpleString(std::string("data['a']=np.matrix("   mat_d.asString(32,'[',']',',')   ",dtype=np.complex128)").c_str());
//pass the inverted matrix to Python
PyRun_SimpleString(std::string("data['b_c']=np.matrix("   mat_d_i.asString(32,'[',']',',')   ",dtype=np.complex128)").c_str());
//inverse in numpy
PyRun_SimpleString(std::string("data['b_p']=np.linalg.inv(data['a'])").c_str());
//flatten the matrices to make comparing them easier (make them 1-dimensional)
PyRun_SimpleString("data['fb_p']=((data['b_p']).flatten().tolist())[0]");
PyRun_SimpleString("data['fb_c']=((data['b_c']).flatten().tolist())[0]");

//make the comparison. The function compare_floats(f1,f2,t) calls numpy.isclose(f1,f2,rtol=t)
//prec is an integer that takes its value from a template function, where I choose the precision I want based on type
PyRun_SimpleString(std::string("res=list(set([compare_floats(data['fb_p'][i],data['fb_c'][i],1e-"  std::to_string(prec)  ") for i in range(len(data['fb_p']))]))[0]").c_str());
//the set above eliminates repeated True and False. If all results are True, we expect that res=[True], otherwise, the test failed somewhere
PyRun_SimpleString(std::string("res = ((len(res) == 1) and res[0])").c_str());
//Now if res is True, then success
 

Комментарии в коде пошагово описывают процедуру.