# Constrained optimization of a cost function for a matrix

I encountered this problem trying to implement an image processing algorithm.

I'm trying to minimize a cost function for a matrix subject to some constraints. The cost function is as follows:

$$J(U)=\sum_{\mathbf{r}}\Bigl(U(\mathbf{r})-\sum_{\mathbf{s}\in N(\mathbf{r})}w_{\mathbf{rs}} U(\mathbf{s})\Bigr)^2$$

Where $w_{\mathbf{rs}}$ is a weighting function:

$$w_{\mathbf{rs}}\propto\left(1+{1\over \sigma_{\mathbf{r}}^2}(U(\mathbf{r})-\mu_{\mathbf{r}})(U(\mathbf{s})-\mu_\mathbf{r})\right)$$

Each point in the matrix has a value of $J(U)$ associated with it. In the above formulas, $\mathbf{r}$ represents a point and $\mathbf{s}$ represents a neighboring (3 x 3 square) point in the matrix. The cost function for each point is a function of the point and surrounding points. $\mu$ and $\sigma$ are constants for each point that I've already computed.

$U(\mathbf{r})$ is the just the value of a parameter $U$ at that point.

Some points are constrained to fixed values of $U$. The objective is to find values of $U$ for every point that minimize the cost function.

The package I'm working with has SciPy, so I have access to these optimization algorithms if any of them would be helpful. I'm relatively new to this kind of math.

I think the problem can be converted into a sparse system of linear equations but I'm not sure how to do so.

Approximate solutions are fine.

• I added parentheses to clarify what $w_{rs}$ is proportional to - if this isn't what you meant, feel free to revert edit. – Zubin Mukerjee Jun 2 '17 at 14:00
• @ZubinMukerjee thanks, that's correct. – AdrianDole Jun 2 '17 at 14:01
• Are there any inequality-type constraints? It looks like the constraints are just that some of the points in $U$ must have specific values. In that case, this reduces to optimization of a multivariable function $J(x_1, x_2, ...)$ with the $x_i$ being any entries of $U$ that aren't fixed. This has no inequality constraint - which means we want to find the points where the gradient of your function is equal to the zero vector. If you have $n$ entries in the matrix $U$ whose values are variable, then you will have a system of $n$ equations after setting the gradient of $J$ equal to $0$. – Zubin Mukerjee Jun 2 '17 at 14:24
• This video has a great visual explanation of using the gradient to minimize/maximize multivariable functions. Because there are no inequality constraints, it shouldn't be necessary to use any other methods such as Karush-Kuhn-Tucker or Lagrance multipliers. Anyway, if there are no inequality-type constraints then I'll flesh out the above comment into an answer. – Zubin Mukerjee Jun 2 '17 at 14:25
• It's worth noting that the set of points (for which the gradient of $J$ is equal to the zero vector) include local maxima and minima, and multidimensional equivalents of saddle points. More tests need to be done before you can be sure that a point is a global minimum. – Zubin Mukerjee Jun 2 '17 at 14:29

## 1 Answer

It looks like none of your constraints in this problem are inequalities. This means we can write $J$ as a function of all of the entries of the matrix $U$ whose values are variable (you mentioned that some entries of $U$ are fixed and some are variable).

Suppose there are $n$ entries of $U$ whose values are variable - call them $x_1, x_2, ..., x_n$. Then, using the two equations from your question, you will be able to write out the cost function $J$, as a function of $x_1, x_2, ..., x_n$.

Once you have done this, you will want to find the set of candidate points (points with $n$ coordinates, in this case) for the global minimum of $J$. To do this, you will want to set the gradient of $J$ equal to the zero vector. The gradient of $J$ is itself a vector with $n$ entries, so setting the gradient of $J$ equal to the zero vector will result in a system of $n$ equations.

The solution set to this system of $n$ equations will be the candidate points for the global minimum of $J$. Some of these points may be local minima, local maxima, or the equivalent of saddle points. We wish to eliminate those points so that only the global minimum remains. You will need to find the all of the second-order partial derivatives of $J$ to compute the Hessian matrix of $J$.

You can then use the method described in the last page of this document to determine which of your candidate points will be your final answer.

The most tedious part of solving this problem will probably actually be writing $J$ in terms of $x_1, x_2, ..., x_n$. It will not be very difficult, but it will be time-consuming.

• Thanks, I will work on implementing this. The function will be working with > 500,000 variables since it's processing images. – AdrianDole Jun 2 '17 at 14:59
• @eliminate1337 Good luck! You will definitely want to look into algorithms for computing Hessian matrices. It would be awful to have a program individually compute $500000^2$ second-order partial derivatives of a $500000$-variable function. – Zubin Mukerjee Jun 2 '17 at 15:02
• @eliminate1337 Out of curiosity, what are you going to use your finished algorithm for? Could you just use existing image-processing code? People have probably already optimized those algorithms to a great extent – Zubin Mukerjee Jun 2 '17 at 15:06
• It's for colorizing an image. It's a relatively new and novel algorithm so there's not much code out there for it. The authors of the paper I'm referencing provided code, but it's a brute force solution. I'm hoping to write a better version that uses SciPy's already optimized minimization functions. – AdrianDole Jun 2 '17 at 15:25
• Very interesting - good luck :) – Zubin Mukerjee Jun 2 '17 at 18:27