How to assess repeatability of multivariate and method-specific outcomes? Method "A" describes biological samples using multivariate "fingerprints" that consist of about 30 different variables. Different variables show different typical distribution and many of them closely correlate one with another. From prior experience it is assumed that we cannot  transform many of the variables to normal distribution.
Method "B" is designed to be an improved version of method "A" and we wish to compare the repeatabilty of these two methods. If we were dealing with single variable, we would perform independent analyses of several samples and use ANOVA in order to compare within-method  to between-methods variability. But here we are dealing with multivariate outputs and we do not wish to perform one analysis per variable. What are the correct approaches to this question?
Resolution
The answer by gui11aume's  answer, provides useful and valuable information. I will adapt the "downstream application" from gui11aume's  answer following by 7 one-way analyses as suggested by AdamO.
 A: This reminds me of cancer diagnostics, where old gene expression signatures are replaced by newer ones, that are of course supposed to be better. But how to show that they are better?
Here are a couple of suggestions to compare the repeatability of the methods.
1. Use co-inertia analysis (CIA).
CIA should be more advertised, unfortunately it is not widely used (no Wikipedia page for example). CIA is a two-table method that works on the same principle as canonical analysis (CA), which is to look for a pair of linear scores with maximum correlation between two sets of multi-densional measurements. Its advantage over CA is that you can do it even if you have more dimensions than observations. You could measure both methods on the same samples to get two coupled tables of 30 columns and $n$ observations. The first pair of principal components should be strongly correlated (if methods really measure the same thing). If method B is better, the residual variance should be smaller than the residual variance of method A. With this approach you address both the agreement of the methods, and their disagreement, which you interpret as noise.
2. Use a distance.
You could use the Euclidean distance in 30 dimensions between the test and the retest to measure the repeatability of a method. You generate a sample of that score for each method and you can compare the samples with the Wilcoxon test.
3. Use downstream application.
You are probably getting these fingerprints to take a decision, or classify patients or biological material. You can count the agreements vs disagreements between tests and retests for both methods and compare them with the Wilcoxon test.
Method 3 is the simplest, but also the most down to earth. Even for high dimensional inputs, decisions are usually quite simple. And however complex our problem is, bear in mind that statistics is the science of decision.
Regarding the question in your comment.

What about using a robust dimensionality reduction method to reduce the multivariate data to a single dimension and analyzing it?

Dimensionality reduction, however robust, will be associated with a loss of variance. If there is a way to transform your multivariate fingerprint into a single score capturing almost all of its variance, then sure, this is by far the best thing to do. But then why is the fingerprint multivariate in the first place?
I assumed from the context of the OP that the fingerprint is multivariate precisely because it is hard to reduce its dimensionality further without losing information. In that case, their repeatability on a single score does not have to be a good proxy for the overall repeatability, because you may neglect the majority of the variance (close to 29/30 in the worst case).
A: I Assume from your question and comment that the 30 output variables can not (easily) or should not be transformed to a single variate. 
One idea to deal with data of $\mathbf{X_A}^{(n \times p_A)} \leftrightarrow   \mathbf{X_B}^{(n \times p_B)}$ is that you could do regression of $\mathbf{X_A}^{(n \times p_A)} \mapsto   \mathbf{X_B}^{(n \times p_B)}$ and vice versa. Additional knowledge (e.g. that variate $i$ in set A corresponds to variate $i$ also in set B) can help to restrict the mapping model and/or with the interpretation.
So what about multi block PCA (or -PLS) which take this idea further? For these methods, both multivariate fingerprints for the same samples (or same individuals) are analyzed together as independent variables, with or without a third dependent block.
R. Brereton: "Chemometrics for Pattern Recognition" discusses some techniques in the last chapter ("Comparing Different Patterns") and googling will lead you to a number of papers, also introductions. Note that your situations sounds similar to problems where e.g. spectroscopic and genetic measurements are analysed together (two matrices with a row-wise correspondence as opposed to analyzing e.g. time series of spectra where a data cube is analysed).
Here's a paper dealing with multi-block analysis: Sahar Hassani: Analysis of -omics data: Graphical interpretation- and validation tools in multi-block methods.
Also, maybe this is a good starting point into another direction: Hoefsloot et.al., Multiset Data Analysis: ANOVA Simultaneous Component Analysis and Related Methods, in: Comprehensive Chemometrics — Chemical and Biochemical Data Analysis(I don't have access to it, just saw the abstract)
A: 30 one way analyses is certainly an option and would be an ideal "table 2" type of analysis, in which an overall performance is summarized in a logical way. It may be the case that Method B produces the first 20 factors with slightly improved precision whereas the last 10 are wildly more variable. You have the issue of inference using a partially ordered space: certainly if all 30 factors are more precise in B, then B is a better method. But there is "grey" area and with the large number of factors, it's almost guaranteed to show up in practice.
If the objective of this research is to land on a single analysis, it's important to consider the weight of each outcome and their endpoint application. If these 30 variables are used in classification, prediction, and/or clustering of observational data, then I would like to see validation of these results and a comparison of A / B in classification (using something like risk stratification tables or mean percent bias), prediction (using the MSE), and clustering (using something like cross validation). This is the proper way of handling the grey area in which you can't say B is better analytically, but works much better in practice.
A: I will try a multivariate ANOVA based on permutation (PERMANOVA) tests appoach. An ordination analisis (based on result on gradient length analysis) could also help.
A: If you could assume multivariate normality (which you said you could not) you could do a Hotelling T2 test of equality of mean vectors to see if you could claim differences between distributions or not.  However although you can't do that you can still theoretically compare the distributions to see if they differ much.  Divide the 30 dimensional space into rectangular grids.  Use these as 30 dimensional bins. Count the number of vectors falling into each bin and apply a chi square test to see if the distributions look the same. The problem with this suggestion is that it requires judiciously selecting the bins in order to cover the data points in an appropriate way.  Also the curse of dimensionality makes it difficult to identify differences between the multivariate distributions without having a very large number of point in each group. I think suggestions that gui11aume gave are sensible.  I don't think the others are. Since comparing the distributions is not feasible in 30 dimensions with a typical sample some form of valid comparison of the mean vectors would seem to me to be appropriate.
