Сжатая матричная функция для поиска пар

Question

Для набора наблюдений:

[a1,a2,a3,a4,a5]

их попарные расстояния

d=[[0,a12,a13,a14,a15]
   [a21,0,a23,a24,a25]
   [a31,a32,0,a34,a35]
   [a41,a42,a43,0,a45]
   [a51,a52,a53,a54,0]]

Даны в сжатой матричной форме (верхний треугольник из вышеприведенного, рассчитано из scipy.spatial.distance.pdist):

c=[a12,a13,a14,a15,a23,a24,a25,a34,a35,a45]

Вопрос в том, что если у меня есть индекс в сжатой матрице, есть ли функция (предпочтительно в python) f , чтобы быстро определить, какие два наблюдения использовались для их расчета?

f(c,0)=(1,2)
f(c,5)=(2,4)
f(c,9)=(4,5)
...

Я пробовал некоторые решения, но ни одно из них не стоит упоминать: (

Answer 1

Формула для индекса сжатой матрицы:

index = d*(d-1)/2 - (d-i)*(d-i-1)/2 + j - i - 1

Где i - индекс строки, j - индекс столбца, а d - длина строки оригинала.(d X d) верхняя треугольная матрица.

Рассмотрим случай, когда индекс ссылается на крайнюю левую ненулевую запись некоторой строки в исходной матрице.Для всех левых индексов

j == i + 1

, поэтому

index = d*(d-1)/2 - (d-i)*(d-i-1)/2 + i + 1 - i - 1
index = d*(d-1)/2 - (d-i)*(d-i-1)/2

С некоторой алгеброй мы можем переписать это как

i**2 + (1 - 2d)*i + 2*index == 0

Тогда мы можем использовать квадратную формулучтобы найти корни уравнения, и мы будем только заботиться о положительном корне.

Если этот индекс соответствует крайней левой ненулевой ячейке, то мы получим положительное целое число в качестве решения, которое соответствует номеру строки.Тогда найти номер столбца просто арифметически.

j = index - d*(d-1)/2 + (d-i)*(d-i-1)/2 + i + 1

Если индекс не соответствует крайней левой ненулевой ячейке, то мы не найдем целочисленный корень, но мы можем взять словоположительный корень в качестве номера строки.

def row_col_from_condensed_index(d,i):
    b = 1 -2*d 
    x = math.floor((-b - math.sqrt(b**2 - 8*i))/2)
    y = i + x*(b + x + 2)/2 + 1
    return (x,y)  

Если вы не знаете d, вы можете вычислить его по длине сжатой матрицы.

((d-1)*d)/2 == len(condensed_matrix)
d = (1 + math.sqrt(1 + 8*len(condensed_matrix)))/2 

Answer 2

Вы можете найти triu_indices полезным.Например,

In []: ti= triu_indices(5, 1)
In []: r, c= ti[0][5], ti[1][5]
In []: r, c
Out[]: (1, 3)

Обратите внимание, что индексы начинаются с 0. Вы можете настроить его, как вам нравится, например:

In []: def f(n, c):
   ..:     n= ceil(sqrt(2* n))
   ..:     ti= triu_indices(n, 1)
   ..:     return ti[0][c]+ 1, ti[1][c]+ 1
   ..:
In []: f(len(c), 5)
Out[]: (2, 4)

Answer 3

Cleary, функция f, которую вы ищете, требует второго аргумента: размерность матрицы - в вашем случае: 5

Первая попытка:

def f(dim,i): 
  d = dim-1 ; s = d
  while i<s: 
    s+=d ; d-=1
  return (dim-d, i-s+d)

Answer 4

Чтобы завершить список ответов на этот вопрос: быстрая, векторизованная версия ответа fgreggs (как предложил Дэвид Маркс) может выглядеть так:

def vec_row_col(d,i):                                                                
    i = np.array(i)                                                                 
    b = 1 - 2 * d                                                                   
    x = (np.floor((-b - np.sqrt(b**2 - 8*i))/2).astype(int)                                      
    y = (i + x*(b + x + 2)/2 + 1).astype(int)                                                    
    if i.shape:                                                                     
        return zip(x,y)                                                             
    else:                                                                           
        return (x,y) 

Мне нужно было сделать эти вычисления для огромных массивов, и ускорение по сравнению с версией без векторизации (https://stackoverflow.com/a/14839010/3631440) (как обычно) довольно впечатляюще (при использовании IPython% timeit):

import numpy as np
from scipy.spatial import distance

test = np.random.rand(1000,1000)
condense = distance.pdist(test)
sample = np.random.randint(0,len(condense), 1000)

%timeit res = vec_row_col(1000, sample)
10000 loops, best of 3: 156 µs per loop

res = []
%timeit for i in sample: res.append(row_col_from_condensed_index(1000, i))
100 loops, best of 3: 5.87 ms per loop

Это примерно в 1009 * 37 раз быстрее в этом примере!

Answer 5

Для повышения эффективности с помощью numpy.triu_indices
используйте это:

def PdistIndices(n,I):
    '''idx = {} indices for pdist results'''
    idx = numpy.array(numpy.triu_indices(n,1)).T[I]
    return idx

Так что I представляет собой массив индексов.

Однако лучшее решение - реализовать оптимизированный поиск методом грубой силы, скажем, в Fortran:

function PdistIndices(n,indices,m) result(IJ)
    !IJ = {} indices for pdist[python] selected results[indices]
    implicit none
    integer:: i,j,m,n,k,w,indices(0:m-1),IJ(0:m-1,2)
    logical:: finished
    k = 0; w = 0; finished = .false.
    do i=0,n-2
        do j=i+1,n-1
            if (k==indices(w)) then
                IJ(w,:) = [i,j]
                w = w+1
                if (w==m) then
                    finished = .true.
                    exit
                endif
            endif
            k = k+1
        enddo
        if (finished) then
            exit
        endif
    enddo
end function

, затем скомпилировать с помощью F2PY и наслаждатьсянепревзойденная производительность.;)

Answer 6

Вот еще одно решение:

import numpy as np

def f(c,n):
    tt = np.zeros_like(c)
    tt[n] = 1
    return tuple(np.nonzero(squareform(tt))[0])

Answer 7

Это дополнение к ответу phynfo и вашему комментарию.Мне не кажется чистым дизайном выводить размер матрицы из длины сжатой матрицы.Тем не менее, вот как вы можете его вычислить:

from math import sqrt, ceil

for i in range(1,10):
   thelen = (i * (i+1)) / 2
   thedim = sqrt(2*thelen + ceil(sqrt(2*thelen)))
   print "compressed array of length %d has dimension %d" % (thelen, thedim)

Аргумент к внешнему квадратному корню всегда должен быть квадратным целым числом, но sqrt возвращает число с плавающей запятой, поэтому при использовании этого параметра необходимо соблюдать осторожность..