1

I am working with some data for which I am interested in calculating some physical parameters. I have a system of linear equations, which I can write in matrix form as:

$$ \textbf{A} \textbf{x} = \textbf{b} $$

where $\textbf{A}$ is a square matrix containing the coefficients of the linear equations, $\textbf{x}$ is a column vector of length $n$ containing the unknown parameters that I want to calculate, and $\textbf{b}$ is a known vector of length $n$.

Given no other constraints, the solution is to simply calculate $\textbf{x}$ by inverting $\textbf{A}$. However, I have hard inequality constraints on $\textbf{x}$ because of physical reasons of the data in question:

$$ 0 \leq x_i \leq c_i \ \forall \ i=1,2,...n $$

where $x_i$ are the values in $\textbf{x}$ and all $c_i$ are known.

Now, these extra inequality constraints make the problem overdetermined. However, I have the choice to remove data (e.g., remove rows from $\textbf{A}$) because I know a priori that some data are more reliable. Thus, I am hoping that I can make the unconstrained problem underdetermined, but given the hard inequality constraints, make the constrained problem exactly determined and calculate a single unique solution of $\textbf{x}$.

To summarize, my question is: how can I solve a determined system of linear equations subject to inequality constraints on the unknown parameters?

I looked into a few potential techniques like linear programming and bounded variable least squares. However, the goal of those methods is to maximize or minimize some objective function, whereas I want an exact solution to my equations. My gut feeling is that a solution should exist but I don't have the linear algebra background to find it, so I appreciate any help!

These are details on my specific problem that might help with a solution:

  • All values in $\textbf{A}$ are between 0 and 1.
  • The value of the $c_i$ is at most 1.
  • $n=36$
llorente
  • 70

2 Answers2

1

I figured out the solution to my question, so I am posting it here in case others are interested.

The main problem is how to apply bounds on the unknown parameters $\textbf{x}$. The key is to substitute the $x_i$ using a transformation that accepts all real numbers and maps the input onto the bounded range of interest (e.g., $0 \leq x_i \leq c_i$). One way to do this is with a sigmoid function. Sigmoid functions have two horizontal asymptotes, which you can scale linearly to match any lower and upper bound. There are a variety of sigmoid functions you can choose from. Using $\arctan$ as an example:

$$ x_i=c_i \left( \frac{1}{\pi} \arctan{y_i} +\frac{1}{2} \right) $$

You can see in the above equation that $y_i$ can be any real number, but $x_i$ is now bounded between 0 and $c_i$.

This transformation gives us a new problem because now we have a system of nonlinear equations. However, we can solve this using Newton's method to numerically iterate for the $y_i$ on the system of equations (see this post for an example). After solving for the $y_i$ we can apply our transformation to convert them to $x_i$, which will automatically be bounded in the specified range.

I tested this $\arctan$ method with synthetic data and it works well. There is probably an art to choosing an appropriate transformation because some sigmoid functions saturate faster as you move to more extreme arguments.

llorente
  • 70
1

You can use Linear Programming with a dummy objective. The usual approach is to use an objective with all zero coefficients: $$\begin{aligned} \min_x \> & 0^Tx \\ & Ax=b \\ & 0 \le x_i \le u_i \end{aligned}$$

This will find a feasible solution and then stop.

Of course, if you are only interested in feasibility, we can actually use any random objective. If I am reading things correctly, you have a unique solution.

Erwin Kalvelagen
  • 4,367
  • 2
  • 13
  • 16