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\$
3
13
\$\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\$
3
  • 1
    \$\begingroup\$ As the only submission that doesn't use builtins, this should be the highest-voted answer. \$\endgroup\$ Dec 15 '17 at 5:46
  • \$\begingroup\$ +1 for calculating something-related-to-determinant in time less than n!! \$\endgroup\$
    – DELETE_ME
    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\$
2
  • \$\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\$
4
  • \$\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\$
    – DELETE_ME
    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\$ Nov 27 '17 at 16:25
  • 2
    \$\begingroup\$ @cairdcoinheringaahing We should really start putting xnor's proposal into practice. \$\endgroup\$ 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
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