How to use Mata’s optimize() in a loop to do iterative estimation?

[Eivind Bernhardsen]

Starting with a pairwise comparison matrix A=[a_ij] _i,j=1…n, where a_ij>=0, a_ii = 1, and a_ij =1/a_ji, where some elements a_ij are missing, it has been suggested that optimal completion of the missing entries can be found by minimizing λ^max(x) with respect to the missing entries x=(x₁,x₂,…,x_k), where λ^max is the largest eigenvalue of A(x). The solution is found by series of iterations min_xr(λ^max(x_r)) where all missing entries are held at fixed values except for x_{r r=1…k}, whithin each series. In the first series all missing variables are set to 1. The following code successfully finds the minimum λ^max for each step and is illustrated for the missing entry x₁ =A[2,3] after all other missing entries are replaced by 1 in A:

void eb(todo, p, lmax, g, H)
{
A=st_matrix("A")
A[2,3]=p[1] +ln(p[1])-ln(p[1])
A[3,2]=1/p[1]
ev=eigenvalues(A)
lmax=abs(ev[1,1])
}
S = optimize_init()
optimize_init_which(S,"min")
optimize_init_evaluator(S, &eb())
optimize_init_params(S, 1)
p = optimize(S)

However, I need to do this for all k missing entries in A, and then again in series until each of the x'es converge. This means I need to run a loop witin a loop. I don’t seem to get any loop working in Mata, and when I try to loop this in Stata where Mata is called to do the estimation, there are several challenges. One is how to let A[i,j] loop over all missing enties, considering that Mata does not use macros. Another is the fact that the do-file in Stata will end whenever there is a change from the Mata mode back to Stata mode, due to the command "end" in Mata which is read "exit" in the do-file. Also, I need to be able to update matrix A, and (ideally also) the initial parameter with each single estimation. Below I post an example matrix convenient for writing the code. Does anybody have a suggestion on how I should proceed?

* Example generated by -dataex-. To install: ssc install dataex
clear
input float(var1 var2 var3) byte(var4 var5) float(var6 var7 var8)
  1   5   3 7 6   6 .33 .25
 .2   1   . 5 .   3   . .14
.33   .   1 . 3   .   6   .
.14  .2   . 1 . .25   . .13
.17   . .33 . 1   .  .2   .
.17 .33   . 4 .   1   . .17
  3   . .17 . 5   .   1   .
  4   7   . 8 .   6   .   1
end

[German Rodriguez]

Why not optimize with respect to all the missing entries at the same time? Here's one way to do it. You first run the above code to read the data into Stata and then run this Mata code:

clear mata
mata:
    M = st_data(.,.)

    // find missing entries
    ij =  J(0, 2, .)
    for(i = 1; i < rows(M); i++) {
        for(j = i + 1; j <= cols(M); j++) {
                if(missing(M[i,j])) ij = ij \ (i,j)
        }
    }
    ij
    // fill missing entries
    real matrix fill(real vector p) {
        external M, ij
        A = M
        for(k = 1; k <= rows(ij); k++) {
            A[ij[k,1], ij[k,2]] = exp(p[k])
            A[ij[k,2], ij[k,1]] = exp(-p[k])
        }
        return(A)
    }
        
    // criterion function
    void eb(todo, p, lmax, g, H) {
        A = fill(p)
        ev = eigenvalues(A)
        lmax = abs(ev[1,1])
    }

    // optimize
    b = J(1, rows(ij), 0)
    S = optimize_init()
    optimize_init_which(S, "min")
    optimize_init_evaluator(S, &eb())
    optimize_init_params(S, b)
    b = optimize(S)
    round(fill(b), 0.0001)
end

The code (1) copies the data into a Mata matrix, (2) loops to find the missing entries above the diagonal; (if i,j is missing so is j,i), (3) defines a function to fill the missing entries from a vector (in this case with 12 elements) working in the log scale, as the elements must all be positive, (4) the criterion function is then quite simple: it fills the matrix and computes the eigenvalue. Here are my results starting from all ones:

: b = optimize(S)
Iteration 0:   f(p) =  10.896069  
Iteration 1:   f(p) =  9.3508042  
Iteration 2:   f(p) =  9.3065907  
Iteration 3:   f(p) =  9.3063978  
Iteration 4:   f(p) =  9.3063978  

: 
: round(fill(b), 0.0001)
            1        2        3        4        5        6        7        8
    +-------------------------------------------------------------------------+
  1 |       1        5        3        7        6        6      .33      .25  |
  2 |      .2        1    .3294        5   1.7178        3    .4653      .14  |
  3 |     .33   3.0355        1   9.8844        3   4.8383        6    .5752  |
  4 |     .14       .2    .1012        1    .5275      .25    .1429      .13  |
  5 |     .17    .5821      .33   1.8957        1    .9279       .2    .1103  |
  6 |     .17      .33    .2067        4   1.0777        1    .2919      .17  |
  7 |       3   2.1489      .17   6.9976        5   3.4252        1    .4072  |
  8 |       4        7   1.7387        8   9.0659        6   2.4559        1  |
    +-------------------------------------------------------------------------+

[Eivind Bernhardsen]

You are completely right, and the code is very neat too. I just learned a lot about how to use Mata. Thank you for helping me out!