how to treat different kind of errors when using numeric solvers?
I would like some help from experienced users of numeric solvers (like optimize.minimize()). I have a system of vector equations that I want to solve numerically. [cid:6d69c188-653b-4b2b-9f36-69bed2a70a15] How should I combine the errors, so the solver can use the feedback best? The resulting errors of the first two equations are about 10 orders of magnitude smaller then the third. What is a "natural" mathematical way to combine such errors? this is my code: def detect_cable(data: np.ndarray): b_measured = np.average(data, axis=1) sensor_coordinates = np.array([[.265, 0, .382], [0, .712, .764], [0, .712, 0], [.752, .712, .382]]) sensor_num = b_measured.shape[0] # print(b_measured) guess_b_earth = np.array( [1.7921, 4.6676, -1.013]) # https://ngdc.noaa.gov/geomag/calculators/magcalc.shtml#igrfwmm # Latitude: 54° 08' 47" N # Longitude: 8° 48' 50" E guess_n = np.array([0, 0, 1]) guess_a = np.array([.2, 0, 0]) guess_r = np.array([0, 2.7, 0]) guess_i = 1000 p0 = np.hstack([guess_b_earth, guess_n, guess_r, guess_a, guess_i]) def error_func_geometry_vs_measurement(x, b_measured): b_earth = x[0:3] n = x[3:6] r = x[6:9] a = x[9:12] i = x[12] k = i * MU_0 / np.pi for sensor_cnt in range(0, sensor_num): r_relativ = 0 if sensor_cnt: # only do this for sensor 1,2 and 3, not 0. # that means that the algorithm will find the projection of 0's point on to the cable on its own. # The other r_relativ vectors depend on r r_relativ = r - n * np.matmul((r - sensor_coordinates[sensor_cnt]), n) / np.linalg.norm(n) else: # for sensor 0 we take r as the projection point on to the cable r_relativ = r - sensor_coordinates[sensor_cnt] r_magnitude = np.linalg.norm(r_relativ, keepdims=True) # print("r: {}".format(r_relativ)) error_0 = -k * np.cross(a - 2 * (r_relativ.T * np.matmul(r_relativ, a)).T, n) / r_magnitude ** 4 + b_earth - b_measured error_0 = np.linalg.norm(error_0, axis=1) error_0 = np.dot(error_0, error_0) # error_1 = np.dot(n, a) ** 2 # error_2 = np.linalg.norm(np.matmul(r, n)) ** 2 error_3 = error_0 # + error_1 + error_2 # print("error_3 {}, error_0 {}, error_1 {}, error_2 {}".format(error_3, error_0, error_1, error_2)) return error_3
Sorry, i cut of part of the function. def detect_cable(data: np.ndarray): b_measured = np.average(data, axis=1) sensor_coordinates = np.array([[.265, 0, .382], [0, .712, .764], [0, .712, 0], [.752, .712, .382]]) sensor_num = b_measured.shape[0] # print(b_measured) guess_b_earth = np.array( [1.7921, 4.6676, -1.013]) # https://ngdc.noaa.gov/geomag/calculators/magcalc.shtml#igrfwmm # Latitude: 54° 08' 47" N # Longitude: 8° 48' 50" E guess_n = np.array([0, 0, 1]) guess_a = np.array([.2, 0, 0]) guess_r = np.array([0, 2.7, 0]) guess_i = 1000 p0 = np.hstack([guess_b_earth, guess_n, guess_r, guess_a, guess_i]) def error_func_geometry_vs_measurement(x, b_measured): b_earth = x[0:3] n = x[3:6] r = x[6:9] a = x[9:12] i = x[12] k = i * MU_0 / np.pi for sensor_cnt in range(0, sensor_num): r_relativ = 0 if sensor_cnt: # only do this for sensor 1,2 and 3, not 0. # that means that the algorithm will find the projection of 0's point on to the cable on its own. # The other r_relativ vectors depend on r r_relativ = r - n * np.matmul((r - sensor_coordinates[sensor_cnt]), n) / np.linalg.norm(n) else: # for sensor 0 we take r as the projection point on to the cable r_relativ = r - sensor_coordinates[sensor_cnt] r_magnitude = np.linalg.norm(r_relativ, keepdims=True) # print("r: {}".format(r_relativ)) error_0 = -k * np.cross(a - 2 * (r_relativ.T * np.matmul(r_relativ, a)).T, n) / r_magnitude ** 4 + b_earth - b_measured error_0 = np.linalg.norm(error_0, axis=1) error_0 = np.dot(error_0, error_0) # error_1 = np.dot(n, a) ** 2 # error_2 = np.linalg.norm(np.matmul(r, n)) ** 2 error_3 = error_0 # + error_1 + error_2 # print("error_3 {}, error_0 {}, error_1 {}, error_2 {}".format(error_3, error_0, error_1, error_2)) return error_3 res = minimize(error_func_geometry_vs_measurement, p0, args=b_measured) x = res.x b_earth = x[0:3] n = x[3:6] r = x[6:9] a = x[9:12] i = x[12] print("B_earth: {}, norm: {}".format(b_earth, np.linalg.norm(b_earth))) print("r: {}, Depth: {}, distance {}".format(r, r[1]-sensor_coordinates[1,1], np.linalg.norm(r))) print("a {}, norm: {}, dot(a,n) {}".format(a, np.linalg.norm(a), np.dot(a,n))) print("n {}, norm: {}, dot(r,n) {}".format(n, np.linalg.norm(a), np.dot(r,n))) print("current I: {}".format(i)) ________________________________ Von: Schuldei, Andreas Gesendet: Donnerstag, 11. November 2021 14:49:07 An: scipy-user@python.org Betreff: how to treat different kind of errors when using numeric solvers? I would like some help from experienced users of numeric solvers (like optimize.minimize()). I have a system of vector equations that I want to solve numerically. [cid:6d69c188-653b-4b2b-9f36-69bed2a70a15] How should I combine the errors, so the solver can use the feedback best? The resulting errors of the first two equations are about 10 orders of magnitude smaller then the third. What is a "natural" mathematical way to combine such errors? this is my code: def detect_cable(data: np.ndarray): b_measured = np.average(data, axis=1) sensor_coordinates = np.array([[.265, 0, .382], [0, .712, .764], [0, .712, 0], [.752, .712, .382]]) sensor_num = b_measured.shape[0] # print(b_measured) guess_b_earth = np.array( [1.7921, 4.6676, -1.013]) # https://ngdc.noaa.gov/geomag/calculators/magcalc.shtml#igrfwmm # Latitude: 54° 08' 47" N # Longitude: 8° 48' 50" E guess_n = np.array([0, 0, 1]) guess_a = np.array([.2, 0, 0]) guess_r = np.array([0, 2.7, 0]) guess_i = 1000 p0 = np.hstack([guess_b_earth, guess_n, guess_r, guess_a, guess_i]) def error_func_geometry_vs_measurement(x, b_measured): b_earth = x[0:3] n = x[3:6] r = x[6:9] a = x[9:12] i = x[12] k = i * MU_0 / np.pi for sensor_cnt in range(0, sensor_num): r_relativ = 0 if sensor_cnt: # only do this for sensor 1,2 and 3, not 0. # that means that the algorithm will find the projection of 0's point on to the cable on its own. # The other r_relativ vectors depend on r r_relativ = r - n * np.matmul((r - sensor_coordinates[sensor_cnt]), n) / np.linalg.norm(n) else: # for sensor 0 we take r as the projection point on to the cable r_relativ = r - sensor_coordinates[sensor_cnt] r_magnitude = np.linalg.norm(r_relativ, keepdims=True) # print("r: {}".format(r_relativ)) error_0 = -k * np.cross(a - 2 * (r_relativ.T * np.matmul(r_relativ, a)).T, n) / r_magnitude ** 4 + b_earth - b_measured error_0 = np.linalg.norm(error_0, axis=1) error_0 = np.dot(error_0, error_0) # error_1 = np.dot(n, a) ** 2 # error_2 = np.linalg.norm(np.matmul(r, n)) ** 2 error_3 = error_0 # + error_1 + error_2 # print("error_3 {}, error_0 {}, error_1 {}, error_2 {}".format(error_3, error_0, error_1, error_2)) return error_3
On Thu, Nov 11, 2021 at 9:07 AM Schuldei, Andreas < andreas.schuldei@th-luebeck.de> wrote:
I would like some help from experienced users of numeric solvers (like optimize.minimize()).
I have a system of vector equations that I want to solve numerically.
How should I combine the errors, so the solver can use the feedback best? The resulting errors of the first two equations are about 10 orders of magnitude smaller then the third. What is a "natural" mathematical way to combine such errors?
The first two equations seem to be constraints (which ought to be satisfied exactly, up to numerical error), while the last one is a measurement with some amount of uncertainty. Notionally, you don't want to just jam them all together into a combined optimization problem where you add up all of those errors together. You want to have a constrained optimization problem where you try to minimize the error of the measurement on the condition that the constraints are satisfied. Right now, you are using a "natural" parameterization where the domain of optimization allows each of these vectors to freely take values all over RR^3 even though they are constrained to a lower-dimensional manifold embedded in the larger space. My usual move here is to try to reparameterize the domain such that the constraints are automatically satisfied (up to numerical error) and just minimize the measurement error. `n` is constrained to be perpendicular to the sensor coordinate vector so instead of 3 free parameters, you really only have one, the azimuthal angle around the `R0` vector (and it only shows up as a unit vector, so there is no magnitude). `a` is constrained to be perpendicular to `n`, so it has 2 degrees of freedom conditioned on the value of `n` (because it does have a magnitude). The optimization domain after reparameterization is not always easy to optimize over. For some kinds of constraints, like these spherical geometry problems, there are specialized optimizers that help optimize over the manifold directly (e.g. Pymanopt). There are also general-purpose constrained optimizers in `scipy.optimize` that let you keep the "natural" parameterization and just express the constraint. These equality constraints are more difficult than inequality constraints, but the `SLSQP` and `trust-constr` minimizes in scipy do allow them. -- Robert Kern
I solved the problem by reducing the dimensions of my input vector to/from the solver and incorporating the constraints into the error functions. The search does not yet find a solution, I assume it is still slightly buggy. But I feel more confident that I am on the right track now. I do like my solution, so I could leave it here for posterity to profit, as well. The trick is to express the right-angled vectors in terms of the first vector and a vector the second angel is rotated around the first one. On StackOverflow, I found a continuous, numerically stable way to find a perpendicular vector to the first one, as well as a nice implementation of the Rodrigues rotation formula. Those two combined should do what I need. ________________________________ Von: Robert Kern <robert.kern@gmail.com> Gesendet: Donnerstag, 11. November 2021 19:37:51 An: SciPy Users List Betreff: [SciPy-User] Re: how to treat different kind of errors when using numeric solvers? On Thu, Nov 11, 2021 at 9:07 AM Schuldei, Andreas <andreas.schuldei@th-luebeck.de<mailto:andreas.schuldei@th-luebeck.de>> wrote: I would like some help from experienced users of numeric solvers (like optimize.minimize()). I have a system of vector equations that I want to solve numerically. [cid:17d10348a1ff456b1e51] How should I combine the errors, so the solver can use the feedback best? The resulting errors of the first two equations are about 10 orders of magnitude smaller then the third. What is a "natural" mathematical way to combine such errors? The first two equations seem to be constraints (which ought to be satisfied exactly, up to numerical error), while the last one is a measurement with some amount of uncertainty. Notionally, you don't want to just jam them all together into a combined optimization problem where you add up all of those errors together. You want to have a constrained optimization problem where you try to minimize the error of the measurement on the condition that the constraints are satisfied. Right now, you are using a "natural" parameterization where the domain of optimization allows each of these vectors to freely take values all over RR^3 even though they are constrained to a lower-dimensional manifold embedded in the larger space. My usual move here is to try to reparameterize the domain such that the constraints are automatically satisfied (up to numerical error) and just minimize the measurement error. `n` is constrained to be perpendicular to the sensor coordinate vector so instead of 3 free parameters, you really only have one, the azimuthal angle around the `R0` vector (and it only shows up as a unit vector, so there is no magnitude). `a` is constrained to be perpendicular to `n`, so it has 2 degrees of freedom conditioned on the value of `n` (because it does have a magnitude). The optimization domain after reparameterization is not always easy to optimize over. For some kinds of constraints, like these spherical geometry problems, there are specialized optimizers that help optimize over the manifold directly (e.g. Pymanopt). There are also general-purpose constrained optimizers in `scipy.optimize` that let you keep the "natural" parameterization and just express the constraint. These equality constraints are more difficult than inequality constraints, but the `SLSQP` and `trust-constr` minimizes in scipy do allow them. -- Robert Kern
participants (2)
-
Robert Kern -
Schuldei, Andreas