Measuring 'parallelness' of vectors I wish to construct some sort of 'measure' (not in the formal sense) of the 'parallelness' of a finite set of $m$ vectors $S = \{v_1 , \ldots , v_m\} \subset \mathbb{R}^n$. This parallelness $p$ should have the following properties:


*

*If $v \in \mathbb{R}^n$ and $\{\lambda_1 , \ldots , \lambda_m \} \subset \mathbb{R}^+$ then
$$
p(\lambda_1 v , \ldots , \lambda_m v) = 1
$$
since, with our set of positive $\lambda$'s, all vectors of the form $\lambda_iv$ point in the 'same diection'

*In any other case
$$
p < 1
$$
to indicate that these vectors are not totally parallel


An easy way to construct such a thing for $m=2$ is using the dot product. Denoting unit vectors by a hat then
$$
p(v_1,v_2) = \hat{v}_1 \cdot \hat{v}_2
$$
Note that $p(v_1,v_2) \leq 1$ because $\hat{v}_1$ and $\hat{v}_2$ are unit vectors.
For more vectors it gets trickier. I currently have the approach
$$
p(v_1, \ldots , v_m) = \left\lVert \frac {\hat{v}_1 + \ldots + \hat{v}_m}{m} \right\rVert
$$
which has the nice property that if the vectors are uniformly distributed over a sphere then $p=0$, no parallelness.
This version is inspired by the mean of circular quantities
Is there a general approach and theory behind what I'm trying to do? Is there a 'better' way to measure how parallel a set of vectors are?
Motivation:
This question is inspired by numerical computations, where I get a number of vector fields and I need to know if the vector fields are parallel. Of course there will be some error in the computation, and so I need to check if the vector fields exceed some level of parallelness.
However, I am interested in whether there is some general approach to get something like a 'standard deviation of direction' or similar in high dimensional space, an analytical tool to approach this kind of problem.
Edit (13/Nov/2017): After considering Raskolnikov's answer, it turns out that I want to first determine whether the vectors are 'parallel' without caring whether they are aligned or anti-aligned, so at this stage $v$ and $-v$ are considered the same, this gives me a 'region' in my vector field. I then want to identify the type of region by comparing all vectors to the first and determining whether they are parallel or anti-parallel. This second step is trivial, it is the first step I am addressing in this question. I therefore update my required properties to be:


*

*If $v \in \mathbb{R}^n$ and $\{\lambda_1 , \ldots , \lambda_m \} \subset \mathbb{R}$ then
$$
p(\lambda_1 v , \ldots , \lambda_m v) = 1
$$
since all vectors of the form $\lambda_iv$ are parallel/antiparallel

*In any other case
$$
p < 1
$$
to indicate that these vectors are not totally parallel/antiparallel

 A: Here's an approach inspired by dbx' comment. We look at the sum of the squared scalar products of a normed vector $\hat{u}$ with the normed vectors $\hat{v}_i$ and we try to maximize this over the possible $\hat{u}$.
$$\max_{\hat{u}}\sum_{i=1}^m (\hat{u} \cdot \hat{v}_i)^2$$
Why choose this measure? Because it has the nice property that it represents a sum of squares of cosines of the angles. So for angles $0$ and $\pi$, a single term is maximal, hence it really checks parallelism and not orientation. 
The other reason is that this variational problem can be translated into an eigenvalue problem. Indeed, we can notate the scalar product as follows
$$\hat{u} \cdot \hat{v}_i = \mathbf{u}^{\text{T}}\mathbf{v}_i$$
where $\mathbf{u}$ and $\mathbf{v}_i$ are column vectors containing our components. Then our optimization problem becomes
$$\max_{\hat{u}}\sum_{i=1}^m (\mathbf{u}^{\text{T}}\mathbf{v}_i)^2 = \max_{\hat{u}}\sum_{i=1}^m (\mathbf{u}^{\text{T}}\mathbf{v}_i)(\mathbf{v}_i^{\text{T}}\mathbf{u})$$
where in the last step, I used the fact that $\mathbf{u}^{\text{T}}\mathbf{v}_i=\mathbf{v}_i^{\text{T}}\mathbf{u}$. Further rearranging we obtain
$$\max_{\hat{u}}\sum_{i=1}^m \mathbf{u}^{\text{T}}(\mathbf{v}_i\mathbf{v}_i^{\text{T}})\mathbf{u} = \max_{\hat{u}} \mathbf{u}^{\text{T}}(\sum_{i=1}^m \mathbf{v}_i\mathbf{v}_i^{\text{T}})\mathbf{u} $$
This optimization problem can be shown to be equivalent with looking for the eigenvectors and eigenvalues of the matrix $\sum_{i=1}^m \mathbf{v}_i\mathbf{v}_i^{\text{T}}$, more particularly the largest eigenvalue is the solution of the problem. (You can work this out for instance through Lagrange optimization, don't forget to put the condition that $\mathbf{u}^{\text{T}}\mathbf{u}=1$.)
$$\sum_{i=1}^m \mathbf{v}_i\mathbf{v}_i^{\text{T}}\mathbf{u}=\lambda \mathbf{u}$$
Interestingly, you can also multiply the equation on the left by $\mathbf{v}_j^T$ and get
$$\sum_{i=1}^m (\mathbf{v}_j^T\mathbf{v}_i)(\mathbf{v}_i^{\text{T}}\mathbf{u})=\lambda (\mathbf{v}_j^T\mathbf{u})$$
which allows you to rephrase the eigenvalue problem into the eigenvalue problem of another matrix $G$ such that $G_{ji}=\mathbf{v}_j^T\mathbf{v}_i=\hat{v}_j\cdot\hat{v}_i$. This is known as the Gram or Gramian matrix and is well-studied in the literature. It contains all the scalar products of the vectors. 
If all $\mathbf{v}_i$ are parallel, the Gram matrix will be all $1$'s or $-1$'s. But in such a way that the rows are all multiples of one another with a factor $\pm 1$. The eigenvalues are thus all zero except one, and that one should be equal to the trace of the Gram matrix, $m$. Thus the result of our optimization will obviously be $m$. So if you want your measure to be $1$, divide by $m$. If all $\mathbf{v}_i$ are mutually orthogonal, the Gram matrix is just the identity matrix, the largest eigenvalue being $1$. Thus your measure will be $1/m$. 
Another interesting property of the Gram matrix: its determinant measures the square of the volume of the parallelotope spanned by the vectors $\hat{v}_i$. If this volume is zero, some of the vectors are parallel (not necessarily all).  
