Eigenvector Concern in Large Fractal-Like Matrix (Possible Bug?)

This is my first post and I have quite a long question. Thanks to any who read.

To start with, I have an operator m=Transpose[n] on a 7 dimensional vector space V where n is constructed as follows:

v = {1, 0, 0, 0, a, 0, b, a, c, 0, g, d, e, g};
n = NestList[Delete[Prepend[#, 0], -1] &, v, 13];

I desire to calculate a basis of the eigenspace for the induced map of m on the 7th exterior power V. To this end I created the induced matrix M3 with the following code:

h[{i_, j_, k_, o_, p_, u_, y_}] :=n[[i]]\[TensorWedge]n[[j]]\[TensorWedge]n[[k]]\[TensorWedge]n[[
o]]\[TensorWedge]n[[p]]\[TensorWedge]n[[u]]\[TensorWedge]n[[y]];

OrderedTuples = Flatten[Table[{i, j, k, o, p, u, y}, {i, 8}, {j, i + 1, 9}, {k,
j + 1, 10}, {o, k + 1, 11}, {p, o + 1, 12}, {u, p + 1, 13}, {y,
u + 1, 14}], 6];

M = Map[h, OrderedTuples];

M2 = Map[Drop[SymmetrizedArrayRules[#], -1] &, M];

Rules = {};

For[i = 1, i < Length[M] + 1, i++, Print[i];
ComposeRule[rulelist_] := rulelist /. HoldPattern[Rule[x_,y_]] :> (x /.Thread[OrderedTuples -> Thread[{i, Range[Length[M]]}]]) -> y;
AppendTo[Rules, Map[ComposeRule, M2[[i]]]]];

M3 = Transpose[Normal[SparseArray[Flatten[Rules], {Length[OrderedTuples], Length[OrderedTuples]}, 0]]];

It takes about an hour to make M3 on my machine (very rough estimate).

Truly, I am interested in the intersection of the eigenspaces (The only eigenvalue of M3 is 1 and if you plot the entries of M3 you see it is like a Sierpenski Gasket, hence the fractal part of the title) over all values of the parameters a,b,c,d,e,g. To start off chose some fixed values like

a = -1236;
b = 5423;
c = 2134;
d = -5453;
e = 2234;
f = 4325;
g = -237;

because sufficiently general parameters should catch all the eigenvectors. To calculate the eignespace I used:

Eigen = DeleteDuplicates[Eigenvectors[N[M3]]];

and obtained a list of 14 eigenvectors, all of which are basis elements. To see this observe that

Length[Eigen]

14

For[i = 1, i < Length[Eigen] + 1, i++,
Print[DeleteDuplicates[Eigen[[i]]]]]

{0.,1.}

{0.,1.}

{0.}

{0.,1.}

{0.,1.}

{0.,1.}

{0.,1.}

{0.,1.}

{0.,1.}

{0.,1.}

{0.,1.}

{0.,1.}

{0.,1.}

{0.,1.}

and

For[i = 1, i < Length[Eigen] + 1, i++, Print[Count[Eigen[[i]], 1.]]]

1

1

0

1

1

1

1

1

1

1

1

1

1

1

Also, I can see exactly which basis vectors these are with

Pos = Sort[Flatten[Map[Position[#, 1.] &, Eigen]]]

{3092, 3311, 3388, 3394, 3396, 3417, 3421, 3423, 3424, 3429, 3430, 3431, 3432}

Now here is the interesting part: From the context of the problem (Algebraic Geometry) I know that

E = SparseArray[{{1603, 1} -> 1, {1606, 1} -> 1, {1695, 1} -> -1, {1698, 1} -> -1}, {3432, 1}];

Should be an eigenvector of M3 and indeed,

ArrayRules[M3.Normal[E]]

{{1603, 1} -> 1, {1606, 1} -> 1, {1695, 1} -> -1, {1698, 1} -> -1, {_, _} -> 0}

as I thought it should be. Is this a bug? Why isn't mathematica producing a basis for the eigenspace? Also, I exported the matrix to matlab and it also calculated the eigenvectors wrong and gave essentially the same answer.

Thanks again!

Ian

Edit 1: I wanted to point out that when I do the same procedure for much smaller cases mathematica is able to give me results that make sense in the bigger theory and this is why it appears to be a scaling issue. Also I have added the following picture to show the fractal structure Edit 2: I have used

Export["M3.mtx", SparseArray[M3]]

as J.M. requested. To create a compressed sparse array version of M3. Here is the link: https://www.dropbox.com/s/8fg8rmtrv494fve/M3.mtx?dl=0

• Ian, you put together a nice question and thorough background, but I wonder if you might be able to extract a minimal example that reproduces the issue you are experiencing. As it stands, the sheer length and complexity of your problem might decrease your chances to get help. – MarcoB Nov 9 '15 at 19:30
• Hi MarcoB, thanks for the feedback. I have been trying to construct a similar example using rotation matrices just to see if the scale is the problem. So far this is to no avail. – Ian Nov 9 '15 at 19:50
• If Mathematica and MATLAB both have the same problem, it is possibly an MKL issue. Which version of both packages did you test? – Oleksandr R. Nov 9 '15 at 20:23
• Mathematica 10.2 uses MKL 11.2 Update 2. Do you know which version MATLAB is using? Also, it would be helpful if you could state which processor you have. Bug fixes for MKL 11.2 are listed here, but nothing seems relevant to this issue. Please try another version of Mathematica or MATLAB, if possible, or alternatively Compress your matrix M3 and upload it somewhere for others to work with, rather than waiting 1 hour to generate it? – Oleksandr R. Nov 9 '15 at 20:31
• Sorry I was wrong about my last post I will edit when I can. I use 10.0.1 for mathematica and 8.6 for MATLAB. My processor is 2.5 GHz Intel Core i5. – Ian Nov 9 '15 at 20:34

You are new, and this is a difficult problem you have been working on. So I'll try to point out some ways to make things easier.

(1) Try to make the code as simple as possible. Little things help. I'll show a few in a moment.

(2) Don't use things like n for matrices.

(3) Don't use upper cases things that look like built-in functions, for example Rules, as variables.

(4) For that matter, don't use upper cases things that look like built-in functions, for example ComposeRule, as a function.

(5) Maybe consider modestly descriptive names for variables, e.g. vec and mat. Not essential, but certainly helpful when you want debugging advice from others.

(6) Try to avoid "global" variables in functions you are defining. Again, not really essential but helpful for others trying to debug.

The next few lines set up the matrix and TensorWedge definitions.

vec = {1, 0, 0, 0, a, 0, b, a, c, 0, g, d, e, g};
mat = NestList[Prepend[Most[#], 0] &, vec, 13];
hh[ilist : {_Integer ..}, m_] := TensorWedge @@ m[[ilist]]

After considerable wrangling with your definitions, attempting on smaller examples, I think I decided that you are taking all 7-minors of the 14×14 original matrix, and placing them in your new matrix. That is to say, I believe your M3 is going to be the same as dwarves (aka "seven miners") below. (Among other things, this means we don't even need the definition for hh above, or any other explicit use of TensorWedge. All the churning through sparse and structured array innards likewise disappears.)

Timing[dwarves = Minors[mat, 7];]

(* Out= {29.161219, Null} *)

Dimensions[dwarves]

(* Out= {3432, 3432} *)

Now for eigenstructure. Notice I made one value approximate so everything will be done in machine arithmetic.

rep = {a -> -1236, b -> 5423., c -> 2134, d -> -5453, e -> 2234,
f -> 4325, g -> -237};

Timing[es = Eigensystem[dwarves /. rep];]

(* Out= {31.953408, Null} *)

We did this using machine double libraries. So we might or might not get something viable. What I obtain below appears to have 13 independent eigenvectors.

evals = es[];
Length[evals]
Union[evals]

(* Out= 3432

Out= {1.} *)

evecs = es[];
maxes = Map[Max[Abs[#]] &, evecs];
posns = Position[maxes, Except[0.], {1}, Heads -> False]

(* Out= {{1}, {2}, {3}, {4}, {9}, {10}, {12}, {16}, {37}, {39},
{45}, {122}, {341}} *)

So that's where I am thus far. If I can make more sense of the eigencomputations I'll edit with an update. Not sure anything has gone wrong though, unless general theory suggests 13 would be, well, an unlucky count.

• On the bright side, you can now get the wrong result much faster. – Daniel Lichtblau Nov 10 '15 at 2:15
• Is it bad that I upvoted because of dwarves? :) – J. M. will be back soon Nov 10 '15 at 2:17
• Serious suggestion 1 (which I can't try just yet). Change the replacement rule to be exact. Use nmat=N[...,100] on the matrix. If this does not exhaust your RAM (or it does, but maybe you also have (crudity/vulgarity warning duly posted) a good goat, then try Eigensystem[nmat,20] and see if the higher precision causes more independent vectors to show up. – Daniel Lichtblau Nov 10 '15 at 2:37
• Serious suggestion 2: You know already the only eigenvalue is 1. Say we call the matrix dwarves, with exact values substituted smat.Try to find the exact NullSpace[smat-IdentityMatrix[Length[smat]]]. If that is too much in terms of computational load, I can come up with a skinny SVD approach next. – Daniel Lichtblau Nov 10 '15 at 2:42
• @J.M. I'll ask Grumpy. – Daniel Lichtblau Nov 10 '15 at 2:46