11
\$\begingroup\$

Given a square matrix, output the matrix's eigenvalues. Each eigenvalue should be repeated a number of times equal to its algebraic multiplicity.

The eigenvalues of a matrix A are scalar values λ such that, for some column vector v, A*v = λ*v. They are also the solutions to the characteristic polynomial of A: det(A - λ*I) = 0 (where I is the identity matrix with the same dimensions as A).

Outputs must be accurate to 3 significant digits. All inputs and outputs will be within the representable range of numerical values for your chosen language.

Builtins are acceptable, but you are encouraged to include solutions that do not use builtins.

Test Cases

In these test cases, I represents the imaginary unit. Complex numbers are written in the form a + b*I. All outputs have 3 significant digits of precision.

[[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]
\$\endgroup\$
12
\$\begingroup\$

Haskell, 576 554 532 507 bytes

No built-ins!

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)

Try it online!

Many thanks @ØrjanJohansen for a total of -47 bytes!

Explanation

First this computes the characteristic polynomial with the Faddeev–LeVerrier algorithm which is the function f. Then the function z computes all roots of that polynomial by iterating g which implements Laguerre's Method for finding a root, once a root is found it is removed and g gets called again until the polynomial has degree 1 which is trivially solved by z[a,b]=[-b/a].

Ungolfed

I re-inlined the functions sum,length,magnitude, fromIntegral, zipWith and (&) as well as the little helper (!). The function faddeevLeVerrier corresponds to f, roots to z and g to laguerre respectively.

-- 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
\$\endgroup\$
  • 1
    \$\begingroup\$ As the only submission that doesn't use builtins, this should be the highest-voted answer. \$\endgroup\$ – Esolanging Fruit Dec 15 '17 at 5:46
  • \$\begingroup\$ +1 for calculating something-related-to-determinant in time less than n!! \$\endgroup\$ – user202729 Dec 15 '17 at 7:18
  • \$\begingroup\$ Thanks guys! @user202729: Initially I oversaw the ! and was really confused :D \$\endgroup\$ – ბიმო Dec 15 '17 at 17:34
6
\$\begingroup\$

Octave, 4 bytes

@eig

Try it online!

Only two bytes more than the MATL golfing language equivalent!

Defines an anonymous function handle to the eig built-in. Interestingly, the MATLAB design philosophy goes against many high-end languages, which like to use DescripteFunctionNamesTakingArguments(), whereas MATLAB and consequently Octave tend to get the shortest unambiguous function name possible. For example, to get a subset of eigenvalues (e.g., the smallest n in absolute magnitude), you use eigs.

As a bonus, here's a function (works in MATLAB, and in theory could work in Octave but their solve is not really up to the task) that doesn't use built-ins, but instead symbolically solves the eigenvalue problem, det(A-λI)=0, and converts it to numerical form using vpa

@(A)vpa(solve(det(A-sym('l')*eye(size(A)))))
\$\endgroup\$
3
\$\begingroup\$

MATL, 2 bytes

Yv

Try it online!

Explanation

I followed the usual advice in numerical linear algebra: instead of writing your own function, use a built-in that is specifically designed to avoid numerical instabilities.

Incidentally, it's shorter. ¯\_(ツ)_/¯

\$\endgroup\$
  • \$\begingroup\$ This begs the question, how long would it be without Yv? \$\endgroup\$ – Sanchises Nov 27 '17 at 10:33
  • \$\begingroup\$ @Sanchises I'm not sure. I'd probably go about it by finding the roots (ZQ) of the characteristic polynomial. But explicitly computing the coeffients of the polynomial may be a lot of work \$\endgroup\$ – Luis Mendo Nov 27 '17 at 11:56
2
\$\begingroup\$

Mathematica, 11 bytes

Eigenvalues

Try it online!

\$\endgroup\$
  • \$\begingroup\$ Yes, I expected a built-in answer before clicking "1 new answer to this question". Let's wait for some non-builtin answer... / Basically Mathematica solution is often <first word in the title> \$\endgroup\$ – user202729 Nov 27 '17 at 9:20
  • \$\begingroup\$ The shortest non-pure built-in I got is First@Eigensystem@#& (20 bytes) \$\endgroup\$ – Mr. Xcoder Nov 27 '17 at 13:19
  • 7
    \$\begingroup\$ I actually agree with user202729 here. While it is fun to joke about Mathematica having a builtin for everything, it is very annoying both as a challenge poster and answering to see a builtin answer equivalent to something you tried relatively hard to do. Golfing (IMO) is about attempting to find the shortest algorithm and implementation of said algorithm, but a builtin answer takes that away from the "sport". \$\endgroup\$ – caird coinheringaahing Nov 27 '17 at 16:25
  • 2
    \$\begingroup\$ @cairdcoinheringaahing We should really start putting xnor's proposal into practice. \$\endgroup\$ – Martin Ender Nov 29 '17 at 14:03
1
\$\begingroup\$

R, 22 bytes

function(m)eigen(m)$va

Try it online!

Takes m as a matrix. Frustratingly, the eigen function in R returns an object of class eigen, which has two fields: values, the eigenvalues, and vectors, the eigenvectors.

However, more annoyingly, the optional argument only.values returns a list with two fields, values, containing the eigenvalues, and vectors, set to NULL, but since eigen(m,,T) is also 22 bytes, it's a wash.

\$\endgroup\$
1
\$\begingroup\$

Julia, 12 bytes

n->eig(n)[1]

Try it online!

Unfortunately, eig returns both the eigenvalues and eigenvectors, as a tuple, so we waste another 9 bytes to lambdify it and grab the first item.

\$\endgroup\$
0
\$\begingroup\$

Python + numpy, 33 bytes

from numpy.linalg import*
eigvals
\$\endgroup\$
0
\$\begingroup\$

Pari/GP, 19 bytes

m->mateigen(m,1)[1]

Try it online!

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.