# What is the best tool to use for weighting matrices

**URL:** <https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218>\
**Category:** Optimization (Mathematical)\
**Created:** [September 15, 2021, 5:16pm UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218 "2021-09-15T17:16:11Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nenndul](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nenndul/32/22149_2.png) [@Nenndul](https://discourse.julialang.org/u/Nenndul)\
**Post date:** [September 15, 2021, 5:16pm UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/1 "2021-09-15T17:16:12Z")

</div>

I have a problem that involves combining several matrices with a weighted average to optimize for a given feature, such as the average value of over the data in the resulting matrix.

The setup is essentially like this.

```julia
weights =(1,2,8)

matrixA =[1 2
          5 8]

matrixB =[5 7
          9 0]

matrixC =[4 5
          7 8]

averageMatrix = (matrixA*weights[1] + matrixB*weights[2] + matrixC*weights[3])/3

matrixAveragePre = accumulate(+, averageMatrix[i,j] for i in 1:2, j in 1:2)

matrixAveragePost = matrixAveragePre[4]/4

```

With these dummy values the output is the expected:

> 20.833333333333332

an output which I would like to optimize by adjusting the weights to produce either a specific number such as, say, 30 or to maximize.

This is a fairly straightforward thing to do using Excel’s built-in Solver tool, but I am not sure what the most convenient equivalent here is. I have looked and JuMP and Optim, but they don’t seem to be solving this sort of problem, unless I am mistaken?

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [September 15, 2021, 7:07pm UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/2 "2021-09-15T19:07:24Z")

</div>

This is easily done in both Optim and JuMP, just define a function that takes `weights` as an input and (in your example) returns `(30 - weights*matrices)^2`, and then minimize this function.

---

<div class="post-metadata">

**Author:** ![Nenndul](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nenndul/32/22149_2.png) [@Nenndul](https://discourse.julialang.org/u/Nenndul)\
**Post date:** [September 15, 2021, 10:25pm UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/3 "2021-09-15T22:25:31Z")

</div>

Thanks, that gets me in the door, but now I’ve set it up as follows:

```julia
function f(weights)

    (30 - weights*matrixAveragePost)^2

end

Optim.minimum(f, weights)

```

which yields: 620.8402777777778, which isn’t much use.

`show(weights)` returns the initial [1.0, 2.0, 8.0]

Am I doing something wrong here? I need the adjusted weights returned.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [September 15, 2021, 11:17pm UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/4 "2021-09-15T23:17:18Z")

</div>

This is easy to do in JuMP:

```nohighlight
using JuMP, Ipopt
matrices = [
    [1 2; 5 8], 
    [5 7; 9 0], 
    [4 5; 7 8],
]
model = Model(Ipopt.Optimizer)
@variable(model, 0 <= weights[1:length(matrices)])
@expression(model, X, sum(weights[i] * matrices[i] for i in 1:length(matrices)))
@objective(model, Min, (30 - X[2, 2])^2)
optimize!(model)
value.(weights)
value.(X)

```

---

<div class="post-metadata">

**Author:** ![Nenndul](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nenndul/32/22149_2.png) [@Nenndul](https://discourse.julialang.org/u/Nenndul)\
**Post date:** [September 16, 2021, 2:07am UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/5 "2021-09-16T02:07:49Z")

</div>

Thanks, this seems to be doing the right thing when I run it in my program. I will need to figure the syntax out a bit, especially how you got the formula for @objective, but I’ll scratch my head on it a bit. At the moment I am just using the the length of the actual matrices rather than the 2’s there.

> Number of nonzeros in equality constraint Jacobian…: 0  
> Number of nonzeros in inequality constraint Jacobian.: 0  
> Number of nonzeros in Lagrangian Hessian…: 0
> 
> Total number of variables…: 5  
> variables with only lower bounds: 5  
> variables with lower and upper bounds: 0  
> variables with only upper bounds: 0  
> Total number of equality constraints…: 0  
> Total number of inequality constraints…: 0  
> inequality constraints with only lower bounds: 0  
> inequality constraints with lower and upper bounds: 0  
> inequality constraints with only upper bounds: 0
> 
> iter objective inf\_pr inf\_du lg(mu) ||d|| lg(rg) alpha\_du alpha\_pr ls  
> 0 1.0000000e+00 0.00e+00 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0  
> 1 1.0000000e+00 0.00e+00 1.00e-02 -1.7 2.00e-02 - 9.90e-01 1.00e+00f 1  
> 2 1.0000000e+00 0.00e+00 1.50e-06 -3.8 1.50e-02 - 1.00e+00 1.00e+00f 1  
> 3 1.0000000e+00 0.00e+00 2.51e-14 -8.6 1.67e-03 - 1.00e+00 1.00e+00f 1
> 
> Number of Iterations…: 3
> 
> ```
> (scaled) (unscaled)
> 
> ```
> 
> Objective…: 1.0000000000000000e+00 1.0000000000000000e+00  
> Dual infeasibility…: 2.5059035640133008e-14 2.5059035640133008e-14  
> Constraint violation…: 0.0000000000000000e+00 0.0000000000000000e+00  
> Complementarity…: 1.1704648838917575e-15 1.1704648838917575e-15  
> Overall NLP error…: 2.5059035640133008e-14 2.5059035640133008e-14
> 
> Number of objective function evaluations = 4  
> Number of objective gradient evaluations = 4  
> Number of equality constraint evaluations = 0  
> Number of inequality constraint evaluations = 0  
> Number of equality constraint Jacobian evaluations = 0  
> Number of inequality constraint Jacobian evaluations = 0  
> Number of Lagrangian Hessian evaluations = 1  
> Total CPU secs in IPOPT (w/o function evaluations) = 0.053  
> Total CPU secs in NLP function evaluations = 0.000
> 
> EXIT: Optimal Solution Found.

However, when I `println(weights)`, it gives me:

> VariableRef[weights[1], weights[2], weights[3], weights[4], weights[5]]

The main thing is that I need to go back and actually apply those weights to the input matrices.

In general, if you think these packages are the right ones for the job I’ll have to spend some time figuring them out.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [September 16, 2021, 2:24am UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/6 "2021-09-16T02:24:47Z")

</div>

> `println(weights)`

Try the last two lines of my code.

> especially how you got the formula

It’s the squared difference between what you have and what you want. We want to minimize that.

---

<div class="post-metadata">

**Author:** ![Nenndul](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nenndul/32/22149_2.png) [@Nenndul](https://discourse.julialang.org/u/Nenndul)\
**Post date:** [September 16, 2021, 3:16am UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/7 "2021-09-16T03:16:20Z")

</div>

> Try the last two lines of my code.

Oh, I understand now.

It didn’t occur to me to modify `value.(weights)` to assign it to a variable to use it subsequently. I’m still not used to needing to explicitly calling outputs like that, sorry.

But the result doesn’t seem to be what I’m expecting. I’m getting:

> [1.8750000000249656, 264.0458596788466, 1.875000000025152]

whereas GRG nonlinear in Excel Solver is giving me the correct weights of:

> [7.604777427, 2.606567535, 3.677612558]

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [September 16, 2021, 6:34am UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/8 "2021-09-16T06:34:35Z")

</div>

You can check a proposed solution for weights by using your code in your original post. A function that does this is

```julia
function check_answer(matrixArray, weightArray)
    L = length(matrixArray)
    m, n = size(matrixArray[1])
    averageMatrix = sum(matrixArray[i]*weightArray[i] for i in 1:L)/L

    matrixAveragePre = accumulate(+, averageMatrix[i,j] for i in 1:m, j in 1:n)
    matrixAveragePost = matrixAveragePre[m,n]/(m*n)

    return matrixAveragePost
end

```

If I have

```julia
matrices = [
    [1 2; 5 8], 
    [5 7; 9 0], 
    [4 5; 7 8],
]

```

and `wA = [1,2,8]` I get `check_answer(matrices, wA ) = 20.833`… as expected.  
If I use the answer `wA = [7.604777427, 2.606567535, 3.677612558]`, then I get  
`check_answer(matrices, wA ) = 22.0564`… so not really an answer if wanting to approximate 30.

Using `wA = [1.8750000000249656, 264.0458596788466, 1.875000000025152]` I get 468.330… but this is because Oscar’s answer is solving a slightly different problem. Your problem is more like

```julia
using JuMP, Ipopt
model = Model(Ipopt.Optimizer)

matrices = [
    [1 2; 5 8], 
    [5 7; 9 0], 
    [4 5; 7 8],
]

m, n = size(matrices[1])

@variable(model, 0 <= weights[1:length(matrices)])
@expression(model, X, sum(weights[i] * matrices[i] for i in 1:length(matrices)))
## New: The next line does the same as the "matrixAveragePre = accumulate(+, ..." line 
@expression(model, Y, sum(X[i,j] for i in 1:m for j in 1:n))
## total objective with averaging denominator 3 and 4 respectively dividing Y:
@objective(model, Min, (30 - Y/(3*4))^2)  

optimize!(model)
value.(weights)

```

Typing `Y` in the REPL returns `16 weights[1] + 21 weights[2] + 24 weights[3]`, and a solution of weights this expression approximately equals 30\*3\*4 = 360.

With this in mind, it should be clear that the problem does not have a _unique_ solution, e.g.  
`check_answer(matrices, [1.25906 2.06343 12.3551] ) = 29.999...`  
`check_answer(matrices, [11.0037 3.51295 4.5904] ) = 30.000...`

---

<div class="post-metadata">

**Author:** ![Nenndul](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nenndul/32/22149_2.png) [@Nenndul](https://discourse.julialang.org/u/Nenndul)\
**Post date:** [September 16, 2021, 8:18pm UTC](https://discourse.julialang.org/t/what-is-the-best-tool-to-use-for-weighting-matrices/68218/9 "2021-09-16T20:18:10Z")

</div>

That works!

I think I understand all the individual bits of what you did there well enough to implement it but the syntax is is going to take me a bit longer to unwrap.

The solution doesn’t have to be unique for what I’m doing, but I’ve found that I could herd GRG solver in the right direction by searching for intermediate solutions. That makes it less likely “cheat” to zeroing out some of the weights in my implementation. Perhaps something like that will work here too.

In any event: Thanks to both of you again.