Q:

How do I make this Solve work faster?

I'm quite new to Mathematica and this code is taking forever to work.
Are there easy ways that will make this same code work much more quickly?
HMmatrix[n_] := Module[{A, a, constraints, vars, AAT, equations, sol},
  A = Table[a[i, j], {i, 1, n}, {j, 1, n}];
  constraints = Table[a[i, j]^2 == 1, {i, 1, n}, {j, 1, n}] // Flatten;
  vars = Table[a[i, j], {i, 1, n}, {j, 1, n}] // Flatten;
  AAT = Dot[A, Transpose[A]];
  equations = Join[{AAT == n*IdentityMatrix[n]}, constraints];
  sol = Solve[equations, vars, Integers]; size = Length[sol]; 
  Return[size]]

A:

Easiest way to solve this problem quickly is to very thoroughly re-state the problem. I do not promise that this is a complete solution to the problem, but I hope it provides significant insight into how to approach finding any such complete solution.
From what I understood from the comments, we are looking for matrices consisting of only $-1$ and $+1$ which satisfy the equation:
$M \times M^T = n \times I_n$
Where $M$ is the square matrix in question, $I_n$ is the $n \times n$ identity matrix, and $n$ is the number of rows or columns in $M$.
By the very nature of the problem, the diagonal elements will always turn out to be $n$ after the multiplication for any square matrix consisting of only 1s (plus or minus). Thus, we are concerned only about every non-diagonal element of the result.
Every single one of these elements meets the criteria that the row ($r$) and column ($c$) involved in forming it satisfy $r \cdot c = 0$.
Let us first find all such possible vectors. Assume $n=2$ for our starting case.
n = 2;
t = Tuples[{-1, 1}, {n}];
r = Outer[Dot, t, t, 1];
p = Position[r, 0];

Line by line, we set n = 2, we construct all possible 2-vectors consisting of -1 or 1, we dot product all of these vectors together to form a matrix n, which we then search for every 0 element and store their positions in p.
The matrix we are looking for, $M$, has the property that every row and column is present in t, and every possible multiple of these rows and columns has indices that are present in p together.
That is, if we look at p right now, it consists of:

{{1, 2}, {1, 3}, {2, 1}, {2, 4}, {3, 1}, {3, 4}, {4, 2}, {4, 3}}

If we collect elements which are paired, we will find that matrices formed from them satisfy the original problem. That is to say that t[[{1,2}]], t[[{1,3}]], t[[{2,4}]], and so on are solutions to the n=2 case.
For n=4 it is somewhat more complicated:

{{1, 4}, {1, 6}, {1, 7}, {1, 10}, {1, 11}, {1, 13}, {2, 3}, {2, 
    5}, {2, 8}, {2, 9}, {2, 12}, {2, 14}, {3, 2}, {3, 5}, {3, 8}, {3, 
    9}, {3, 12}, {3, 15}, {4, 1}, {4, 6}, {4, 7}, {4, 10}, {4, 11}, {4, 
    16}, {5, 2}, {5, 3}, {5, 8}, {5, 9}, {5, 14}, {5, 15}, {6, 1}, {6, 
    4}, {6, 7}, {6, 10}, {6, 13}, {6, 16}, {7, 1}, {7, 4}, {7, 6}, {7, 
    11}, {7, 13}, {7, 16}, {8, 2}, {8, 3}, {8, 5}, {8, 12}, {8, 14}, {8,
     15}, {9, 2}, {9, 3}, {9, 5}, {9, 12}, {9, 14}, {9, 15}, {10, 
    1}, {10, 4}, {10, 6}, {10, 11}, {10, 13}, {10, 16}, {11, 1}, {11, 
    4}, {11, 7}, {11, 10}, {11, 13}, {11, 16}, {12, 2}, {12, 3}, {12, 
    8}, {12, 9}, {12, 14}, {12, 15}, {13, 1}, {13, 6}, {13, 7}, {13, 
    10}, {13, 11}, {13, 16}, {14, 2}, {14, 5}, {14, 8}, {14, 9}, {14, 
    12}, {14, 15}, {15, 3}, {15, 5}, {15, 8}, {15, 9}, {15, 12}, {15, 
    14}, {16, 4}, {16, 6}, {16, 7}, {16, 10}, {16, 11}, {16, 13}}

If we just pick randomly from these, we will not get a solution except by chance. However, if we pick in turn so that every new index appears in a pair with every prior chosen index, we can get a solution. One example of such is t[[{1,4,6,7}]]:

{{-1, -1, -1, -1}, {-1, -1, 1, 1}, {-1, 1, -1, 1}, {-1, 1, 1, -1}}

We can automate this process significantly by recognizing this as a graph problem and applying appropriate functions.
g = Graph[Evaluate[Map[DirectedEdge @@ # &, p]]];
c = FindClique[g, {n}, All]

{{10, 11, 13, 16}, {7, 11, 13, 16}, {6, 10, 13, 16}, {6, 7, 13, 
    16}, {4, 10, 11, 16}, {4, 7, 11, 16}, {4, 6, 10, 16}, {4, 6, 7, 
    16}, {9, 12, 14, 15}, {8, 12, 14, 15}, {5, 9, 14, 15}, {5, 8, 14, 
    15}, {3, 9, 12, 15}, {3, 8, 12, 15}, {3, 5, 9, 15}, {3, 5, 8, 
    15}, {2, 9, 12, 14}, {2, 8, 12, 14}, {2, 5, 9, 14}, {2, 5, 8, 
    14}, {2, 3, 9, 12}, {2, 3, 8, 12}, {2, 3, 5, 9}, {2, 3, 5, 8}, {1, 
    10, 11, 13}, {1, 7, 11, 13}, {1, 6, 10, 13}, {1, 6, 7, 13}, {1, 4, 
    10, 11}, {1, 4, 7, 11}, {1, 4, 6, 10}, {1, 4, 6, 7}}

That said, I am not 100% familiar with the mathematics involved here as I do not have much personal experience with graph theory. I am not going to claim that these are all of the solutions to this problem, but it does appear that each of these cliques is an individual solution to the problem.
This should be generally faster than attempting to directly solve larger cases. However, it is still rather slow. The first solution it can find reasonably quickly is: 

{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {-1, -1, -1, -1, 
  -1, -1, 1, 1, 1, 1, 1, 1}, {-1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 
    1}, {-1, -1, 1, -1, 1, 1, -1, 1, 1, -1, -1, 1}, {-1, -1, 1, 1, -1, 
    1, 1, -1, 1, -1, 1, -1}, {-1, -1, 1, 1, 1, -1, 1, 1, -1, 
    1, -1, -1}, {-1, 1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1}, {-1, 1, -1, 
    1, -1, 1, 1, 1, -1, -1, -1, 1}, {-1, 1, -1, 1, 1, -1, -1, 1, 1, -1, 
    1, -1}, {-1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, 1}, {-1, 1, 1, 
    1, -1, -1, -1, -1, 1, 1, -1, 1}, {1, -1, -1, 1, 1, -1, 1, -1, 
    1, -1, -1, 1}}

Interestingly, there are apparently no solutions for $n=6$ or $n=10$. I haven't been able to spend the time to check the $n=14$ case. There are definitely solutions for 2, 4, 8, and 12, though exactly how many there are for 12 is not something I've spent the time to check.

A:

Not a full answer...
I took what proved to be a similar approach as @eyorble, in that I turn it into a graph problem.
I also just demonstrate how to get a lot of solutions, but not the guaranteed total number of solutions. I think the total number blows up combinatorically, so probably not worth counting.
The main idea is that the matrix $A$ is invertible, with vector elements that almost form an orthonormal basis. For $a$ and $b$ in $A$, $a.b$ is either $n$ or $0$.
Let the first vector $a_{1}$ be the vector $[{1,1,....,1}]$ Then every other vector must have an equal number of $1$ and $-1$. Create all of those permutations. Let $n=12$
n = 12; 
pp = Permutations[Flatten@{Table[1, {n/2}], Table[-1, {n/2}]}];

Take the dot product of all of them, but turning it into an adjacency matrix where it is a $1$ if the dot product is zero, and zero otherwise.
gg = Outer[If[Dot[#1, #2] != 0, 0, 1] &, pp, pp, 1];//Timing
(* {0.875, Null} *)

Now find 2000 cliques (meaning 2000 solutions for $A$) that have exactly 11 nodes.
res = FindClique[AdjacencyGraph@gg, {n - 1}, 2000]; // Timing

(* {0.125, Null} *)

Length@res
(* 2000 *)

If you just run FindClique[AdjacencyGraph@gg] it takes much longer, and the largest clique also has 11 elements.
There are no answers with this method for $n=14$, and it blows up on $n=16$. I am pretty sure that there's a constructive method to get to higher dimensions without having to do a brute force search, as in this method.

