[SciPy-User] linalg question: unique sign of qr ?

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

[SciPy-User] linalg question: unique sign of qr ?

josef.pktd
I'm trying to do QR updating and my updating algorithm doesn't manage to match the signs of numpy or scipy linalg.qr.


Is there a sign convention for R matrix of QR?

Wikipedia says " If A is invertible, then the factorization is unique if we require that the diagonal elements of R are positive."

scipy and numpy seems to have all diagonal values negative except for the last one which is positive.

R'R looks to be the same even with different signs (and the same as x dot x of the data). So similar to eigenvalue decomposition we need to pick some signs.

Extra question: If I want to match the signs, do I need to change signs by column or by rows in the upper triangular R matrix?

Josef


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: linalg question: unique sign of qr ?

josef.pktd



On Wed, Sep 10, 2014 at 11:38 AM, <[hidden email]> wrote:
I'm trying to do QR updating and my updating algorithm doesn't manage to match the signs of numpy or scipy linalg.qr.


Is there a sign convention for R matrix of QR?

Wikipedia says " If A is invertible, then the factorization is unique if we require that the diagonal elements of R are positive."

scipy and numpy seems to have all diagonal values negative except for the last one which is positive.

R'R looks to be the same even with different signs (and the same as x dot x of the data). So similar to eigenvalue decomposition we need to pick some signs.

Extra question: If I want to match the signs, do I need to change signs by column or by rows in the upper triangular R matrix?


and for comparison cholesky has positive diagonal elements

>>> np.round(np.linalg.qr(xy, mode='r') / np.linalg.cholesky(xy.T.dot(xy)).T , 3)
array([[ -1.,  -1.,  -1.,  -1.,  -1.],
       [ nan,  -1.,  -1.,  -1.,  -1.],
       [ nan,  nan,  -1.,  -1.,  -1.],
       [ nan,  nan,  nan,  -1.,  -1.],
       [ nan,  nan,  nan,  nan,   1.]])

my updating QR seems to have random signs, that can be fixed by sign adjustments by rows

>>> np.round(res[-1] * np.sign(res[-1].diagonal())[:,None] / np.linalg.cholesky(xy.T.dot(xy)).T , 3)
array([[  1.,   1.,   1.,   1.,   1.],
       [ nan,   1.,   1.,   1.,   1.],
       [ nan,  nan,   1.,   1.,   1.],
       [ nan,  nan,  nan,   1.,   1.],
       [ nan,  nan,  nan,  nan,   1.]])


BTW: updating the R of QR without keeping track of Q after new observations have arrived is just two lines (after finding a reference).

Josef
 

Josef



_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: linalg question: unique sign of qr ?

Geordie McBain-4
2014-09-11 1:57 GMT+10:00  <[hidden email]>:

>
>
>
> On Wed, Sep 10, 2014 at 11:38 AM, <[hidden email]> wrote:
>>
>> I'm trying to do QR updating and my updating algorithm doesn't manage to
>> match the signs of numpy or scipy linalg.qr.
>>
>>
>> Is there a sign convention for R matrix of QR?
>>
>> Wikipedia says " If A is invertible, then the factorization is unique if
>> we require that the diagonal elements of R are positive."
>>
>> scipy and numpy seems to have all diagonal values negative except for the
>> last one which is positive.
>>
>> R'R looks to be the same even with different signs (and the same as x dot
>> x of the data). So similar to eigenvalue decomposition we need to pick some
>> signs.
>>
>> Extra question: If I want to match the signs, do I need to change signs by
>> column or by rows in the upper triangular R matrix?

Hi Josef,
  I found the same nonunicity of numpy.linalg.qr.  My fix was
similarly to enforce the positive diagonal on R.

If A = QR and S = sgn diag R, then S is an involution (S inv S = I)
and QR = Q S inv S R = (Q S) (inv S R), and the last is another QR
decomposition but with the property required for uniqueness.  The code
I use is

  Q, R = np.linalg.qr(A)
  Q = Q * np.sign(np.diag(R))

HTH

G. D. McBain
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: linalg question: unique sign of qr ?

josef.pktd



On Wed, Sep 10, 2014 at 7:50 PM, Geordie McBain <[hidden email]> wrote:
2014-09-11 1:57 GMT+10:00  <[hidden email]>:
>
>
>
> On Wed, Sep 10, 2014 at 11:38 AM, <[hidden email]> wrote:
>>
>> I'm trying to do QR updating and my updating algorithm doesn't manage to
>> match the signs of numpy or scipy linalg.qr.
>>
>>
>> Is there a sign convention for R matrix of QR?
>>
>> Wikipedia says " If A is invertible, then the factorization is unique if
>> we require that the diagonal elements of R are positive."
>>
>> scipy and numpy seems to have all diagonal values negative except for the
>> last one which is positive.
>>
>> R'R looks to be the same even with different signs (and the same as x dot
>> x of the data). So similar to eigenvalue decomposition we need to pick some
>> signs.
>>
>> Extra question: If I want to match the signs, do I need to change signs by
>> column or by rows in the upper triangular R matrix?

Hi Josef,
  I found the same nonunicity of numpy.linalg.qr.  My fix was
similarly to enforce the positive diagonal on R.

If A = QR and S = sgn diag R, then S is an involution (S inv S = I)
and QR = Q S inv S R = (Q S) (inv S R), and the last is another QR
decomposition but with the property required for uniqueness.  The code
I use is

  Q, R = np.linalg.qr(A)
  Q = Q * np.sign(np.diag(R))

Thanks Geordie for the explanation and confirming this.

It would have saved me an hour of hunting for non-existing bugs had I thought about this earlier.

Josef

 

HTH

G. D. McBain
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user