Wartości własne macierzy

11

Biorąc pod uwagę macierz kwadratową, wyprowadzaj wartości własne macierzy. Każda wartość własna powinna być powtarzana kilka razy równa jej algebraicznej krotności.

Te wartości własnych macierzy Asą wartości skalarne λ, tak że przez pewien wektor kolumny v, A*v = λ*v. Są to również rozwiązania charakterystycznego wielomianu A: det(A - λ*I) = 0(gdzie Ijest macierz tożsamości o takich samych wymiarach jak A).

Dane wyjściowe muszą mieć dokładność do 3 cyfr znaczących. Wszystkie dane wejściowe i wyjściowe będą w reprezentatywnym zakresie wartości liczbowych dla wybranego języka.

Wbudowane są dopuszczalne, ale zachęca się do uwzględnienia rozwiązań, które nie używają wbudowanych.

Przypadki testowe

W tych przypadkach testowych Ireprezentuje wyimaginowaną jednostkę. Liczby zespolone są zapisywane w formularzu a + b*I. Wszystkie wyjścia mają 3 znaczące cyfry precyzji.

[[42.0]] -> [42.0]
[[1.0, 0.0], [0.0, 1.0]] -> [1.00, 1.00]
[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]] -> [16.1, -1.12, -1.24e-15]
[[1.2, 3.4, 5.6, 7.8], [6.3, 0.9, -5.4, -2.3], [-12.0, -9.7, 7.3, 5.9], [-2.5, 7.9, 5.3, 4.4]] -> [7.20 + 5.54*I, 7.20 - 5.54*I, -4.35, 3.75]
[[-3.22 - 9.07*I, 0.193 + 9.11*I, 5.59 + 1.33*I, -3.0 - 6.51*I, -3.73 - 6.42*I], [8.49 - 3.46*I, -1.12 + 6.39*I, -8.25 - 0.455*I, 9.37 - 6.43*I, -6.82 + 8.34*I], [-5.26 + 8.07*I, -6.68 + 3.72*I, -3.21 - 5.63*I, 9.31 + 3.86*I, 4.11 - 8.82*I], [-1.24 + 9.04*I, 8.87 - 0.0352*I, 8.35 + 4.5*I, -9.62 - 2.21*I, 1.76 - 5.72*I], [7.0 - 4.79*I, 9.3 - 2.31*I, -2.41 - 7.3*I, -7.77 - 6.85*I, -9.32 + 2.71*I]] -> [5.18 + 16.7*I, -24.9 - 2.01*I, -5.59 - 13.8*I, 0.0438 - 10.6*I, -1.26 + 1.82*I]
[[-30.6 - 73.3*I, 1.03 - 15.6*I, -83.4 + 72.5*I, 24.1 + 69.6*I, 52.3 + 2.68*I, 23.8 + 98.0*I, 96.8 + 49.7*I, -26.2 - 5.87*I, -52.4 + 98.2*I, 78.1 + 6.69*I], [-59.7 - 66.9*I, -26.3 + 65.0*I, 5.71 + 4.75*I, 91.9 + 82.5*I, -94.6 + 51.8*I, 61.7 + 82.3*I, 54.8 - 27.8*I, 45.7 + 59.2*I, -28.3 + 78.1*I, -59.9 - 54.5*I], [-36.0 + 22.9*I, -51.7 + 10.8*I, -46.6 - 88.0*I, -52.8 - 32.0*I, -75.7 - 23.4*I, 96.2 - 71.2*I, -15.3 - 32.7*I, 26.9 + 6.31*I, -59.2 + 25.8*I, -0.836 - 98.3*I], [-65.2 - 90.6*I, 65.6 - 24.1*I, 72.5 + 33.9*I, 1.47 - 93.8*I, -0.143 + 39.0*I, -3.71 - 30.1*I, 60.1 - 42.4*I, 55.6 + 5.65*I, 48.2 - 53.0*I, -3.9 - 33.0*I], [7.04 + 0.0326*I, -12.8 - 50.4*I, 70.1 - 30.3*I, 42.7 - 76.3*I, -3.24 - 64.1*I, 97.3 + 66.8*I, -11.0 + 16.5*I, -40.6 - 90.7*I, 71.5 - 26.2*I, 83.1 - 49.4*I], [-59.5 + 8.08*I, 74.6 + 29.1*I, -65.8 + 26.3*I, -76.7 - 83.2*I, 26.2 + 99.0*I, -54.8 + 33.3*I, 2.79 - 16.6*I, -85.2 - 3.64*I, 98.4 - 12.4*I, -27.6 - 62.3*I], [82.6 - 95.3*I, 55.8 - 73.6*I, -49.9 + 42.1*I, 53.4 + 16.5*I, 80.2 - 43.6*I, -43.3 - 3.9*I, -2.26 - 58.3*I, -19.9 + 98.1*I, 47.2 + 62.4*I, -63.3 - 54.0*I], [-88.7 + 57.7*I, 55.6 + 70.9*I, 84.1 - 52.8*I, 71.3 - 29.8*I, -3.74 - 19.6*I, 29.7 + 1.18*I, -70.6 - 10.5*I, 37.6 + 99.9*I, 87.0 + 19.0*I, -26.1 - 82.0*I], [69.5 - 47.1*I, 11.3 - 59.0*I, -84.3 - 35.1*I, -3.61 - 35.7*I, 88.0 + 88.1*I, -47.5 + 0.956*I, 14.1 + 89.8*I, 51.3 + 0.14*I, -78.5 - 66.5*I, 2.12 - 53.2*I], [0.599 - 71.2*I, 21.7 + 10.8*I, 19.9 - 97.1*I, 20.5 + 37.4*I, 24.7 + 40.6*I, -82.7 - 29.1*I, 77.9 + 12.5*I, 94.1 - 87.4*I, 78.6 - 89.6*I, 82.6 - 69.6*I]] -> [262. - 180.*I, 179. + 117.*I, 10.3 + 214.*I, 102. - 145.*I, -36.5 + 97.7*I, -82.2 + 89.8*I, -241. - 104.*I, -119. - 26.0*I, -140. - 218.*I, -56.0 - 160.*I]
Mego
źródło
Powiązane ? Prawdopodobnie związany ? (zależy od podejścia)
użytkownik202729

Odpowiedzi:

12

Haskell , 576 554 532 507 bajtów

Brak wbudowanych!

import Data.Complex
s=sum
l=length
m=magnitude
i=fromIntegral
(&)=zip
t=zipWith
(x!a)b=x*a+b
a#b=[[s$t(*)x y|y<-foldr(t(:))([]<$b)b]|x<-a]
f a|let c=[1..l a];g(u,d)k|m<-[t(+)a b|(a,b)<-a#u&[[s[d|x==y]|y<-c]|x<-c]]=(m,-s[s[b|(n,b)<-c&a,n==m]|(a,m)<-a#m&c]/i k)=snd<$>scanl g(0<$c<$c,1)c
p?x|let f=foldl1(x!);c=l p-1;n=i c;q p=init$t(*)p$i<$>[c,c-1..];o=f(q p)/f p;a|d<-sqrt$(n-1)*(n*(o^2-f(q$q p)/f p)-o^2)=n/last(o-d:[o+d|m(o-d)<m(o+d)])=last$p?(x-a):[x|m a<1e-9]
z[a,b]=[-b/a]
z p=p?0:z(init$scanl1(p?0!)p)

Wypróbuj online!

Wielkie dzięki @ ØrjanJohansen za łącznie -47 bajtów!

Wyjaśnienie

Najpierw oblicza się charakterystyczny wielomian za pomocą algorytmu Faddeeva – LeVerriera, który jest funkcją f. Następnie funkcja zoblicza wszystkie pierwiastki wielomianu przez iterację, gktóra implementuje metodę Laguerre'a znajdowania pierwiastka, po znalezieniu pierwiastka jest on usuwany i gwywoływany ponownie, aż wielomian ma stopień 1, który jest trywialnie rozwiązany z[a,b]=[-b/a].

Nie golfił

I ponownie inlined funkcji sum, length, magnitude, fromIntegral, zipWithi (&), a także małego pomocnika (!). W przypadku funkcji faddeevLeVerrierodpowiada f, rootsdo zi gz laguerreodpowiednio.

-- Transpose a matrix/list
transpose a = foldr (zipWith(:)) (replicate (length a) []) a

-- Straight forward implementation for matrix-matrix multiplication
(#) :: [[Complex Double]] -> [[Complex Double]] -> [[Complex Double]]
a # b = [[sum $ zipWith (*) x y | y <- transpose b]|x<-a]


-- Faddeev-LeVerrier algorithm
faddeevLeVerrier :: [[Complex Double]] -> [Complex Double]
faddeevLeVerrier a = snd <$> scanl go (zero,1) [1..n]
  where n = length a
        zero = replicate n (replicate n 0)
        trace m = sum [sum [b|(n,b)<-zip [1..n] a,n==m]|(m,a)<-zip [1..n] m]
        diag d = [[sum[d|x==y]|y<-[1..n]]|x<-[1..n]]
        add as bs = [[x+y | (x,y) <- zip a b] | (b,a) <- zip as bs]
        go (u,d) k = (m, -trace (a#m) / fromIntegral k)
          where m = add (diag d) (a#u)


-- Compute roots by succesively removing newly computed roots
roots :: [Complex Double] -> [Complex Double]
roots [a,b] = [-b/a]
roots   p   = root : roots (removeRoot p)
  where root = laguerre p 0
        removeRoot = init . scanl1 (\a b -> root*a + b)

-- Compute a root of a polynomial p with an initial guess x
laguerre :: [Complex Double] -> Complex Double -> Complex Double
laguerre p x = if magnitude a < 1e-9 then x else laguerre p new_x
  where evaluate = foldl1 (\a b -> x*a+b)
        order' = length p - 1
        order  = fromIntegral $ length p - 1
        derivative p = init $ zipWith (*) p $ map fromIntegral [order',order'-1..]
        g  = evaluate (derivative p) / evaluate p
        h  = (g ** 2 - evaluate (derivative (derivative p)) / evaluate p)
        d  = sqrt $ (order-1) * (order*h - g**2)
        ga = g - d
        gb = g + d
        s = if magnitude ga < magnitude gb then gb else ga
        a = order /s
        new_x = x - a
ბიმო
źródło
1
Jako jedyne zgłoszenie, które nie korzysta z wbudowanych poleceń, powinno to być najwyżej ocenione głosowanie.
Esolanging Fruit
+1 za obliczenie czegoś związanego z wyznacznikiem w czasie krótszym niż n!!
user202729,
Dzięki chłopaki! @ user202729: Początkowo nadzorowałem !i byłem bardzo zdezorientowany: D
ბიმო
6

Oktawa , 4 bajty

@eig

Wypróbuj online!

Tylko dwa bajty więcej niż odpowiednik języka golfowego MATL!

Definiuje anonimowy uchwyt funkcji do eigwbudowanego. Co ciekawe, filozofia projektowania MATLAB jest sprzeczna z wieloma wysokiej klasy językami, które lubią używać DescripteFunctionNamesTakingArguments(), podczas gdy MATLAB i w konsekwencji Octave mają tendencję do uzyskiwania możliwie najkrótszej jednoznacznej nazwy funkcji. Na przykład, aby uzyskać ów ubset wartości własnych (np najmniejszy nw absolutnej wielkości), należy użyć eigs.

Jako bonus, oto funkcja (działa w MATLAB, i teoretycznie może działać w Octave, ale solvetak naprawdę nie jest do zadania), która nie używa wbudowanych, ale zamiast tego symbolicznie rozwiązuje problem wartości własnej det(A-λI)=0i przekształca go do postaci numerycznej za pomocąvpa

@(A)vpa(solve(det(A-sym('l')*eye(size(A)))))
Sanchises
źródło
3

MATL , 2 bajty

Yv

Wypróbuj online!

Wyjaśnienie

Postępowałem zgodnie ze zwykłymi wskazówkami z numerycznej algebry liniowej: zamiast pisać własną funkcję, użyj wbudowanej funkcji zaprojektowanej specjalnie w celu uniknięcia niestabilności numerycznej.

Nawiasem mówiąc, jest krótszy. ¯ \ _ (ツ) _ / ¯

Luis Mendo
źródło
To nasuwa pytanie, jak długo by to było bez Yv?
Sanchises
@ Sanchises Nie jestem pewien. Prawdopodobnie zająłbym się tym, znajdując pierwiastki ( ZQ) charakterystycznego wielomianu. Ale wyraźne obliczenie współczynników wielomianu może być bardzo pracochłonne
Luis Mendo
2

Mathematica, 11 bajtów

Eigenvalues

Wypróbuj online!

J42161217
źródło
Tak, oczekiwałem wbudowanej odpowiedzi przed kliknięciem „1 nowa odpowiedź na to pytanie”. Poczekajmy na jakąś niewbudowaną odpowiedź ... / Zasadniczo rozwiązanie Mathematica jest często <pierwszym słowem w tytule>
user202729
Najkrótszy nie-czysty wbudowany, jaki mam First@Eigensystem@#&(20 bajtów)
Pan Xcoder
7
Zgadzam się z użytkownikiem202729 tutaj. Chociaż fajnie jest żartować, że Mathematica ma wbudowane narzędzie do wszystkiego, jest to bardzo denerwujące zarówno jako plakat z wyzwaniem, jak i odpowiedź na pytanie, czy wbudowana odpowiedź jest odpowiednikiem czegoś, co próbowałeś stosunkowo ciężko zrobić. Gra w golfa (IMO) polega na próbie znalezienia najkrótszego algorytmu i implementacji tego algorytmu, ale wbudowana odpowiedź odbiera to „sportowi”.
caird coinheringaahing
2
@cairdcoinheringaahing Naprawdę powinniśmy zacząć stosować propozycję xnor w praktyce .
Martin Ender
1

R , 22 bajty

function(m)eigen(m)$va

Wypróbuj online!

Przyjmuje się mjako matryca. Frustrujące jest to, że eigenfunkcja w R zwraca obiekt klasy eigen, który ma dwa pola: valueswartości własne i vectorswektory własne.

Jednak bardziej denerwujące jest to, że opcjonalny argument only.valueszwraca a listz dwoma polami, valueszawierającymi wartości własne, i vectorsustawiony na NULL, ale ponieważ eigen(m,,T)jest również 22 bajtów, jest to pranie.

Giuseppe
źródło
1

Julia , 12 bajtów

n->eig(n)[1]

Wypróbuj online!

Niestety eigzwraca zarówno wartości własne, jak i wektory własne, jako krotkę, więc marnujemy kolejne 9 bajtów, aby je lambdyfikować i zabrać pierwszy przedmiot.

Uriel
źródło
0

Python + numpy, 33 bajty

from numpy.linalg import*
eigvals
Uriel
źródło