Создать жорданову матрицу из собственных значений, используя NumPy

Question

У меня есть ndarray собственных значений и их кратностей (например, np.array([(2.2, 2), (3, 3), (5, 1)])). Мне нужно вычислить жорданову матрицу для этих собственных значений без использования Python циклов и итераций (списки, for loops и c.), Только с использованием функций NumPy.

Я решил построить матрицу следующим образом:

1. Создать эти блоки, используя np.vectorize и np.eye с np.fill_diagonal:

Объедините блоки в одну матрицу, используя hstack и vstack.

Но у меня есть две проблемы:

1. Вот фрагмент моего кода, создающего код :

def eye(t):
    eye = np.eye(t[1].astype(int),k=1)
    return eye

def jordan_matrix(X: np.ndarray) -> np.ndarray:
    dim = np.sum(X[:,1].astype(int))
    eyes = np.vectorize(eye, signature='(x)->(n,m)')(X)
    return eyes

И я получаю ошибку ValueError: could not broadcast input array from shape (3,3) into shape (2,2)

Мне нужно создать дополнительные нулевые матрицы для заполнения пространства, которое не используется созданными блоками, но их размеры являются переменными, и я не могу понять, как их создать, не используя Python for и его эквиваленты .

Я на правильном пути? Как мне выйти из этой проблемы?

Answer 1

np.vectorize будет в основном l oop под капотами. Мы могли бы использовать NumPy funcs для фактической векторизации на уровне Python. Вот один из таких способов -

def blockwise_jordan(a):
    r = a[:,1].astype(int)
    v = np.repeat(a[:,0],r)
    out = np.diag(v)
    n = out.shape[1]

    fillvals = np.ones(n, dtype=out.dtype)
    fillvals[r[:-1].cumsum()-1] = 0
    out.flat[1::out.shape[1]+1] = fillvals
    return out

Пример выполнения -

In [52]: X = np.array([(2.2, 2), (3, 3), (5, 1)])

In [53]: blockwise_jordan(X)
Out[53]: 
array([[2.2, 1. , 0. , 0. , 0. , 0. ],
       [0. , 2.2, 0. , 0. , 0. , 0. ],
       [0. , 0. , 3. , 1. , 0. , 0. ],
       [0. , 0. , 0. , 3. , 1. , 0. ],
       [0. , 0. , 0. , 0. , 3. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 5. ]])

Оптимизация # 1

Мы можем заменить три последних шага для выполнения условное присвоение единиц и нулей, вот так -

out.flat[1::n+1] = 1
c = r[:-1].cumsum()-1
out[c,c+1] = 0

Answer 2

Вот мое решение:

def jordan(a):
    e = a[:,0]  # eigenvalues
    m = a[:,1].astype('int')  # multiplicities
    d = np.repeat(e, m) # main diagonal
    ones = np.ones(d.size - 1)
    ones[np.cumsum(m)[:-1] -1] = 0
    j = np.diag(d) + np.diag(ones, k=1)
    return j

Редактировать: только что понял, что мое решение почти такое же, как у Divakar.