Extract off-diagonal elements

[Jesse Wursten]

Assume you have a symmetric square matrix sigma

T = 6
sigma = invvech((1 .. T*(T+1)/2)')

Is there an (computationally) efficient way to extract some of the off-diagionals? For example in the example code, you get

I would like to end up with a vector that for example takes the first and second off-diagonal, that is

2 8 13 17 20 3 9 14 18

In my actual purpose, the number can be anything, not necessarily this nice sequence of integers. The matrix can be very large and the operation executed many times, hence why a somewhat efficient solution would be useful. So far, I haven't managed to find one yet, even though I feel it should not be that complicated. Some select(sigma, ???) magic could do it I think, but again, I haven't managed to find the right syntax. E.g. in pseudo code, it could be something along the lines of select if (row-column)<3 & (row-column) > 1 (I think).

[Jesse Wursten]

   
T = 10
    sigma = invvech((1 .. T*(T+1)/2)')
    sigma = sigma[1..T-1, 1..T-1]
    selection = J((T)*(T-1)/2, 1, .)

    count = 1
    for(i=1; i<=T; i++){
        for(j = 1; j<=T-i; j++){
            selection[count++] = j
        }
    }

    selectionCriterion = vec(lowertriangle(invvech(selection))) :>1
    selectionCriterion2 = vec(lowertriangle(invvech(selection))) :<=3
    
    selectionCriteria = selectionCriterion :* selectionCriterion2
    
    sigma
    select(vec(sigma), selectionCriteria)

After plowing on, I've managed to come up with this. It seems to work, but it's not elegant at all and the select is not quite fast. For T = 5000, it took .356s. Constructing the selection criterion only needs to be doen once, but the selection step I might need to do a couple thousand times, which would quickly mean several dozen minutes. Acceptable, but improvement is still desirable.

[Scott Merryman]

This seems to work:

mata:
T = 6
sigma = invvech((1 .. T*(T+1)/2)')
sigma
OffDiag = diagonal(sigma[1..rows(sigma),2..cols(sigma)])  /// 
  \ diagonal(sigma[1..rows(sigma)-2,3..cols(sigma)])
OffDiag
end

[Jesse Wursten]

This seems to work:

mata:
T = 6
sigma = invvech((1 .. T*(T+1)/2)')
sigma
OffDiag = diagonal(sigma[1..rows(sigma),2..cols(sigma)]) ///
\ diagonal(sigma[1..rows(sigma)-2,3..cols(sigma)])
OffDiag
end

Is there a way to scale this such that the number of off-diagonals is a parameter? E.g. depending on user input, I might need 2 of them, or 20.

[Christophe Kolodziejczyk]

Working further on Scott's solution I came up with this

clear mata
mata:

real colvector selOffDiag(real rowvector diag, real matrix X)
{

        real colvector XoffDiag
        real scalar T
        real scalar i
        
        T = rows(X)
        
        XoffDiag = J(0,1,.)
        for (i=1;i <= cols(diag) ; i++)
        {
                XoffDiag = XoffDiag\diagonal(X[|diag[i],1\T,T-diag[i]+1|])
        
        }
        
        return(XoffDiag)

}
T = 6
sigma = invvech((1 :: T*(T+1)/2))
selOffDiag((2,3),sigma)
end
exit

[Christophe Kolodziejczyk]

Otherwise there is this possibility with a loop (works only with symmetric matrices though)

clear mata
mata:


function selOffDiag(real rowvector diag, real matrix X)
{
        
        real scalar K, nEl
        real scalar i, j
        
        
        index = J(0,1,.)
        
        K = rows(X)

        l = 1
        for (j=1 ; j <= cols(diag); j++)
        {
                nEl = K-diag[j]+1
                for (i=1; i<= nEl ; i++)
                {
                    if (i==1) index = index\diag[j]
                    else index = index\index[l-1]+(K-(i-2))
                    l++
                }
        }
        
        return(vech(X)[index])


}

T = 6
sigma = invvech((1 :: T*(T+1)/2))
sigma
selOffDiag((2,4),sigma)

end

[Jesse Wursten]

I'll test it after the weekend to see which works better, thanks so much to all of you for the suggestions!

[Jesse Wursten]

I ended up going for a slightly modified version of #5 as it is twice as fast as my own attempt, thanks to you both!

real colvector extractDiag(real scalar p, real scalar T, real matrix X) {
    real colvector XoffDiag
    
    XoffDiag = J(0,1,.)
    for (i=2; i <= p+1 ; i++)    {
        XoffDiag = XoffDiag\diagonal(X[|i,1\T,T-i+1|])
    }
    
    return(XoffDiag)
}

  1.       8.91 /        1 =     8.913
  2.       15.1 /        1 =    15.125

[Rafal Raciborski (StataCorp)]

Here is a rough draft of a slightly faster Mata function that does not rely on joining matrices. Instead, I preallocate the output vector and fill it in as I go through the input matrix.

The function assumes a square matrix and omits the main diagonal, thus to match Christophe's indices I must substract 1 from mine.

cscript

mata:

transmorphic colvector __offdiag(transmorphic matrix X, real vector L, U) {

    real scalar N, n_L, n_U, NN, ix
    real matrix ixx
    transmorphic vector res
    
    N = rows(X)
    if (cols(X) != N) {
        errprintf("symmetric matrix required\n")
        exit(198)
    }
    
    n_L = length(L)
    n_U = length(U)
    
    if (n_L==1 & L==0) n_L = 0
    if (n_U==1 & U==0) n_U = 0
    
    if (any(L:>=N) | any(U:>=N)) {
        errprintf("diagonal # must be < matrix dimension\n")
        exit(198)
    }
    
    NN = N*n_L + N*n_U - sum(L) - sum(U)
    
    res = J(NN,1,0)
    ix = 0

    // upper off-diagonals
    for (i=n_U; i>=1; i--) {
        ixx = (1::N-U[i]),(U[i]+1::N)
        for (j=1; j<=rows(ixx); j++) res[++ix] = X[ixx[j,1],ixx[j,2]]
    }
    
    // lower off-diagonals
    for (i=1; i<=n_L; i++) {
        ixx = (L[i]+1::N),(1::N-L[i])
        for (j=1; j<=rows(ixx); j++) res[++ix] = X[ixx[j,1],ixx[j,2]]
    }
    return(res)
}

real colvector selOffDiag(real rowvector diag, real matrix X)
{
    real colvector XoffDiag
    real scalar T, i
    T = rows(X)
    XoffDiag = J(0,1,.)
    for (i=1;i <= cols(diag) ; i++) {
        XoffDiag = XoffDiag\diagonal(X[|diag[i],1\T,T-diag[i]+1|])
    }
    return(XoffDiag)
}

end

mata:
    T = 1000
    sigma = round(1000*runiform(T,T))

    timer_clear()
    
    timer_on(1)
    m1 = selOffDiag((2,5),sigma)
    timer_off(1)
    
    timer_on(2)
    m2 = __offdiag(sigma,(1,4),0)
    timer_off(2)
    
    timer()
    
    asserteq(m1,m2)    
end

Timings:

N    selOffDiag()    __offdiag()
------------------------------------
1000          .003        .001
10,000        .184        .005
20,000        .710        .010
30,000       1.607        .015

[Jesse Wursten]

Wow, that is really fast. Is there a way to distribute this such that everyone can use it easily, akin to "ssc install" in non-Mata Stata?

[Jesse Wursten]

Just as a note to anyone finding this thread at a later point - the syntax of __offdiag is __offdiag(<matrix to extract diagonals from>, <vector listing each lower diagonal you want to extract>, <vector listing each upper diagonal you want to extract>).

In particular, if you want off-diagonal 1 up to 4, you should specify (1,2,3,4) as second option. Not (1,4). Of course you can also use (1..4).

[Rafal Raciborski (StataCorp)]

We have just released an update to Stata; the update includes a new undocumented Mata function subdiagget() that extracts subdiagonals from square matrices, see help mata subdiagget().

Note that unlike the __offdiag() function above, subdiagget() extracts the upper diagonals left to right. If you wish to extract the upper diagonals right to left, simply reverse the subscripts:

d = subdiagget(X,0,5..1)

rather than

d = subdiagget(X,0,1..5)