Различия между QR-разложением R и Numpy

Question

Я работаю через большую кодовую базу R (v3.6.0) и пытаюсь понять, что она делает. Для этого я перевожу часть кода R на Python (v3.6.5), используя Numpy (v1.14.3). У меня есть фрагмент кода R, который, кажется, работает нормально:

> v<-c(1,1,1,1)
> qrout<-qr(v)
> qr.Q(qrout)
     [,1]
[1,] -0.5
[2,] -0.5
[3,] -0.5
[4,] -0.5
> qr.R(qrout)
     [,1]
[1,]   -2

Эквивалент Python не отлично:

>>> import numpy as np
>>> v=np.ones(4)
>>> v
array([1., 1., 1., 1.])
>>> np.linalg.qr(v)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/opt/python/3.6.5/lib/python3.6/site-packages/numpy/linalg/linalg.py", line 753, in qr
    _assertRank2(a)
  File "/opt/python/3.6.5/lib/python3.6/site-packages/numpy/linalg/linalg.py", line 195, in _assertRank2
    'two-dimensional' % a.ndim)
numpy.linalg.linalg.LinAlgError: 1-dimensional array given. Array must be two-dimensional

Глядя на документы этоПохоже, что в R используется LAPACK DQRDC(2) / DGEQP3 / ZGEQP3, в то время как Numpy использует LAPACK dgeqrf, zgeqrf, dorgqr и zungqr. Очевидно, что R доволен одномерной матрицей, а Numpy - нет.

ВОПРОС

Как мне повторить QR-факторизацию R, используя Numpy?

Answer 1

Как указано в сообщении об ошибке

Массив должен быть двумерным

In [7]: qr(v[:,None])                                                                     
Out[7]: 
(array([[-0.5],
        [-0.5],
        [-0.5],
        [-0.5]]), array([[-2.]]))

Edit
^{То, что следует, ничем не отличается от приведенного выше кода, но кто знает ...}

In [28]: from numpy.linalg import qr 
    ...: from numpy import ones

In [29]: v = ones(4) ; print(v.shape) ; print(v[:,None].shape) # adding a dimension
(4,)
(4, 1)

In [30]: q, r = qr(v[:, None])

In [31]: print(q) ; print() ; print(r)                 
[[-0.5]
 [-0.5]
 [-0.5]
 [-0.5]]

[[-2.]]

In [32]:

В массивах Python / Numpy может быть только одно измерение, но для qr требуется двумерный массив.

Например, в Python транспозиция не изменяет размеры того, что по сути является одномерным вектором.

In [9]: print(v); print(v.T)                                                              
[1 1 1 1]
[1 1 1 1]

[10]: print(v.shape); print((v.T).shape)                                                
(4,)
(4,)

В R qr() пытается привести свой ввод к двумерному массиву (матрице), поэтому qr() делает этот шаг за вас, в то время как в Python вы должны сделать это явно.

^{Самый идиоматичный способ добавить измерение в массив Numpy - использовать None в объекте среза, чтобы обозначить добавление к нему фиктивного измерения.}