2008-10-28 6 views
21

Je pense que cela doit être simple, mais je ne peux pas faire les choses ...algorithme pour les indices de coefficients matrice triangulaire

J'ai une matrice triangulaire MxM, dont les coefficients sont stockés dans un vecteur, ligne par rangée. Par exemple:

M = [ m00 m01 m02 m03 ] 
     [  m11 m12 m12 ] 
     [   m22 m23 ] 
     [    m33 ] 

est stocké comme

coef[ m00 m01 m02 m03 m11 m12 m13 m22 m23 m33 ] 

Maintenant, je suis à la recherche d'un algorithme non récurrent qui me donne la taille de la matrice M et l'indice de tableau de coefficients i

unsigned int row_index(i,M) 

et

unsigned int column_index(i,M) 

de l'élément de matrice auquel il se réfère. Ainsi, row_index(9,4) == 3, column_index(7,4) == 2 etc. si le comptage d'index est basé sur zéro.

EDIT: Plusieurs réponses utilisant une itération ont été données. Est-ce que quelqu'un connaît des expressions algébriques?

+0

Pouvez-vous reformuler cela pour que les éléments du tableau soient des lettres plutôt que les nombres - il est difficile à 100% de comprendre exactement ce que vos fonctions sont censées faire lorsque vous utilisez à la fois des valeurs de rang/col basées sur zéro et des valeurs basées sur un seul dans le tableau lui-même. – Alnitak

+0

@Alnitak: terminé. –

+0

data [row * num_cols + column] est l'expression classique d'un tableau C 2D stocké en mémoire sous forme de tableau unique. –

Répondre

7

est ici une solution algébrique (la plupart du temps):

unsigned int row_index(unsigned int i, unsigned int M){ 
    double m = M; 
    double row = (-2*m - 1 + sqrt((4*m*(m+1) - 8*(double)i - 7)))/-2; 
    if(row == (double)(int) row) row -= 1; 
    return (unsigned int) row; 
} 


unsigned int column_index(unsigned int i, unsigned int M){ 
    unsigned int row = row_index(i, M); 
    return i - M * row + row*(row+1)/2; 
} 

EDIT: ROW_INDEX fixe()

+0

Je ne pense pas que cela fonctionne pour chaque M. Par exemple, cette solution donne row_index (12,5) = 2, tandis que la bonne réponse (si je l'ai fait correctement) est 3. – mattiast

+0

Mon observation aussi bien. M = 6 donne 3 valeurs erronées pour row_index (et donc aussi pour column_index). –

+0

Ouais, c'était trop facile. Ça fonctionne maintenant. –

0

Il m'a fallu du temps pour comprendre ce dont vous aviez besoin! :)

unsigned int row_index(int i, int m) 
{ 
    int iCurrentRow = 0; 
    int iTotalItems = 0; 
    for(int j = m; j > 0; j--) 
    { 
     iTotalItems += j; 

     if((i+1) <= iTotalItems) 
      return iCurrentRow; 

     iCurrentRow ++; 
    } 

    return -1; // Not checking if "i" can be in a MxM matrix. 
} 

Désolé oublièrent autre fonction .....

unsigned int column_index(int i, int m) 
{ 
    int iTotalItems = 0; 
    for(int j = m; j > 0; j--) 
    { 
     iTotalItems += j; 

     if((i+1) <= iTotalItems) 
      return m - (iTotalItems - i); 
    } 

    return -1; // Not checking if "i" can be in a MxM matrix. 
} 
2

Il doit être que

i == col + row*(M-1)-row*(row-1)/2 

donc une façon de trouver col et la ligne est à itérer valeurs possibles de la ligne:

for(row = 0; row < M; row++){ 
    col = i - row*(M-1)-row*(row-1)/2 
    if (row <= col < M) return (row,column); 
} 

Ceci est au st non récursif, je ne sais pas si vous pouvez le faire sans itération.

Comme on peut le voir à partir de ceci et d'autres réponses, il faut à peu près calculer la ligne pour calculer la colonne, donc il pourrait être intelligent de faire les deux en une seule fonction.

2

Il peut y avoir un intelligent une doublure pour ces derniers, mais (moins toute erreur de vérification):

unsigned int row_index(unsigned int i, unsigned int M){ 
    unsigned int row = 0; 
    unsigned int delta = M - 1; 
    for(unsigned int x = delta; x < i; x += delta--){ 
     row++; 
    } 
    return row; 
} 

unsigned int column_index(unsigned int i, unsigned int M){ 
    unsigned int row = 0; 
    unsigned int delta = M - 1; 
    unsigned int x; 
    for(x = delta; x < i; x += delta--){ 
     row++; 
    } 
    return M + i - x - 1; 
} 
+1

L'index de ligne est OK, l'index de colonne a un décalage de 2. Pourriez-vous changer l'instruction de retour dans 'return M + i - x - i; '? Merci. –

+1

Oups. Je l'ai réparé. –

18

one-liners à la fin de cette réponse, explication suivante :-)

la co tableau efficace a: M éléments pour la première rangée (ligne 0, dans votre indexation), (M-1) pour le deuxième (ligne 1), et ainsi de suite, pour un total de M + (M-1) + & hellip; 1 = M (M + 1)/2 éléments. Il est un peu plus facile de travailler à partir de la fin, car le tableau de coefficients toujours a 1 élément pour la dernière rangée (rangée M-1), 2 éléments pour l'avant-dernière rangée (rangée M-2), 3 éléments pour la rangée M-3, et ainsi de suite.Les dernières rangées K occupent les dernières positions K (K + 1)/2 du tableau de coefficients.

Supposons maintenant que vous avez un index i dans le tableau de coefficients. Il y a M (M + 1)/2-1-i aux positions supérieures à i. Appelez ce numéro i '; vous voulez trouver le plus grand k tel que k (k + 1)/2 ≤ i '. La forme de ce problème est la source même de votre misère - pour autant que je puisse voir, vous ne pouvez pas éviter de prendre des racines carrées :-)

Eh bien faisons-le quand même: k (k + 1) ≤ 2i 'signifie (k + 1/2) - 1/4 ≤ 2i ', ou de manière équivalente k ≤ (sqrt (8i' + 1) -1)/2. Permettez-moi d'appeler le plus grand k tel que K, puis il y a K lignes qui sont plus tard (sur un total de M lignes), de sorte que le row_index (i, M) est M-1-K. Comme pour l'index de la colonne, K (K + 1)/2 éléments de l'i 'sont dans les rangées suivantes, donc il y a j' = i'-K (K + 1)/2 éléments plus tard dans cette rangée (qui a K + 1 éléments au total), donc l'index de la colonne est K-j '. [Ou de manière équivalente, cette ligne commence à (K + 1) (K + 2)/2 à partir de la fin, donc nous avons juste besoin de soustraire la position de départ de cette rangée de i: i- [M (M + 1)/2 - (K + 1) (K + 2)/2]. Il est encourageant de constater que les deux expressions donnent la même réponse]

(pseudo-) Code [en utilisant ii au lieu de i « que certaines langues ne permettent pas de variables avec des noms comme i » ;-)]:

row_index(i, M): 
    ii = M(M+1)/2-1-i 
    K = floor((sqrt(8ii+1)-1)/2) 
    return M-1-K 

column_index(i, M): 
    ii = M(M+1)/2-1-i 
    K = floor((sqrt(8ii+1)-1)/2) 
    return i - M(M+1)/2 + (K+1)(K+2)/2 

Bien sûr, vous pouvez les transformer en interligne en substituant les expressions de ii et K. Vous devrez peut-être faire attention aux erreurs de précision, mais il existe des moyens de trouver la racine carrée entière qui n'exigent pas de calcul en virgule flottante . En outre, pour citer Knuth: "Méfiez-vous des bugs dans le code ci-dessus, je l'ai seulement prouvé correct, pas essayé." Si je peux me permettre de faire une autre remarque: le simple fait de conserver toutes les valeurs dans un tableau M × M prendra au maximum deux fois la mémoire, et en fonction de votre problème, le facteur 2 pourrait être insignifiant par rapport aux améliorations algorithmiques , et pourrait valoir la peine d'échanger le calcul éventuellement coûteux de la racine carrée pour les expressions plus simples que vous aurez. [Edit: BTW, vous pouvez prouver ce plancher ((sqrt (8ii + 1) -1)/2) = (isqrt (8ii + 1) -1)/2 où isqrt (x) = plancher (sqrt (x)) est une racine carrée entière, et la division est une division entière (troncature, la valeur par défaut en C/C++/Java etc.) - donc si vous avez des soucis de précision, vous ne devez vous soucier de l'implémentation d'un nombre entier correct racine carrée.]

+0

C'est exactement ce que j'avais besoin de trouver! Dans mon cas, j'essaie d'effectuer des calculs sur chaque paire d'éléments entre deux vecteurs sur un GPU, donc la ligne et la colonne sont des indices d'éléments. Dans ce cas, il ne s'agit pas seulement de doubler l'empreinte mémoire, même si cela représente également une grande partie du problème. Il y a quelques ajustements que je dois faire pour traiter l'adressage de la mémoire, mais votre explication donne un bon cadre. Merci! – Omegaman

+1

pouvez-vous montrer à quoi ressemblerait l'algèbre pour une matrice triangulaire sans les éléments diagonaux (c'est-à-dire pour les éléments M (M-1)/2)? – frank

+2

@frank Bien sûr, c'est similaire. Les K derniers rangs ont K (K-1)/2 éléments. Donc, étant donné i (indexation basée sur 0), nous voulons le plus grand K tel que K (K-1)/2 ShreevatsaR

1

J'ai réfléchi un peu et j'ai obtenu le résultat suivant. Notez que vous obtenez à la fois une ligne et des colonnes sur un seul coup.

Hypothèses: lignes de départ est un 0. Les colonnes sont à 0. N commencent à 0

Notation

N = taille de la matrice (M était dans le problème initial)

m = indice de la élément

Le code est psuedo

function ind2subTriu(m,N) 
{ 
    d = 0; 
    i = -1; 
    while d < m 
    { 
    i = i + 1 
    d = i*(N-1) - i*(i-1)/2 
    } 
    i0 = i-1; 
    j0 = m - i0*(N-1) + i0*(i0-1)/2 + i0 + 1; 
    return i0,j0 
} 

Et certains octave/matlab code

function [i0 j0]= ind2subTriu(m,N) 
I = 0:N-2; 
d = I*(N-1)-I.*(I-1)/2; 
i0 = I(find (d < m,1,'last')); 
j0 = m - d(i0+1) + i0 + 1; 

Que pensez-vous?

En décembre 2011, un très bon code pour ce faire a été ajouté à GNU/Octave.Potentiellement, ils étendront sub2ind et ind2sub. Le code peut être trouvé pour le moment comme fonctions privées ind2sub_tril et sub2ind_tril

+0

JuanPi: Vous semblez avoir un ancien compte non enregistré, et un nouveau compte enregistré. J'ai marqué le message pour attirer l'attention d'un mod, mais vous devriez poster sur Meta pour demander une fusion de compte. – bdonlan

+0

@JuanPi - Signalez cette question et indiquez vos deux comptes, nous nous ferons un plaisir de les régler. –

2

L'explication de ShreevatsaR est excellente et m'a aidé à résoudre mon problème. Cependant, l'explication et le code fournis pour l'index de colonne donnent une réponse légèrement différente de ce que la question demande.

Pour réitérer, il y a j '= i' - K (K + 1)/2 éléments dans la rangée après i. Mais la rangée, comme tous les autres rangs, a M éléments. L'indice de colonne (basé sur zéro) est donc y = M - 1 - j '.

Le pseudo-code correspondant est:

column_index(i, M): 
    ii = M(M+1)/2-1-i 
    K = floor((sqrt(8ii+1)-1)/2) 
    jj = ii - K(K+1)/2 
    return M-1-jj 

La réponse donnée par ShreevatsaR, K - j », est l'indice de colonne quand on commence à compter (zéro) à partir de la diagonale de la matrice. Par conséquent, son calcul donne column_index (7,4) = 0 et non, comme spécifié dans la question, column_index (7,4) = 2.

+0

pouvez-vous montrer à quoi ressemblerait l'algèbre pour une matrice triangulaire sans les éléments diagonaux (c'est-à-dire pour les éléments M (M-1)/2, également avec des index commençant à 0)? – frank

Questions connexes