@Gauss  Your designation  NullSpace^-1(Y)  does not make sense, since it is neither a matrix or operator, but  a subspace, so just write  NullSpace(Y) .

If  Y  is a known matrix, for example, such as  Y=<< 1 | -2 | 1>, <-6 | 3 | 0>, <4 | -5 | 2 >>, then to find NullSpace(Y) necessary to solve the homogeneous system of linear equations:

restart;

Y:=<<1 | -2 | 1> , <-6 | 3 | 0> , <4 | -5 | 2>>;

LinearAlgebra[LinearSolve](Y, <0,0,0>);

We see that the NullSpace(Y) is one-dimensional subspace of  R^3  spanned by the vector  <1, 2, 3>

@Markiyan Hirnyk  I mean, if  topic-starter wants to obtain a unique solution, but not an infinite set of solutions, expressed in terms of the extra unknowns.

@Markiyan Hirnyk   Yes, of course

convert(exp(x), Sum, x=1, dummy=n);

@Gauss   but for a given matrix  A  the solution is unique, and your  NewMatrix3  is contrary to found  solution.

@PatrickT   I deleted the extra bracket.

Picture of output loaded in another way.

@Gauss   On my machine, the code is executed in a few seconds (in Classic M 12).  Run  restart  command  first.

@toandhsp  No, it can not. AreIsometric  procedure works only with the sets in the Euclidean plane, rather than in space.

@Alejandro Jakubi  Thanks for the explanation of the reason of the error.

@taro yamada  My method is only suitable for Classic Worksheet. If you write

A:=x/2+y/2:

normal(A/a);

then you see    

And then you just replace  symbol  a  by  the empty symbol  `` and get the desired output:

In the previous embodiment, the diagonal entries are always  0 , and should be random integers from  -9  to  9 . New version is free from this drawback:

roll:=rand(-9 .. 9):
Matrix(4, {seq((i, i)=roll(), i=1..4), seq(seq((i, j)=roll()+I*roll(), j=i+1..4), i=1..3)}, shape=hermitian);

Factoring  procedure takes 1/2 sec for this calculation, and the improvements made by Joe Riel and Carl Love even less time. See  http://www.mapleprimes.com/posts/141668-Partitions-Of-A-Natural-Number-Into-Factors

ts:=time():

Factoring(9!):

nops(%);

%%[1..20];

time()-ts;

I have already answered a similar question about the matrix  A . About the matrix  B  all is still easier. There is an inaccuracy in determining  B. Apparently instead of  -1/(2*M*(M-2))   should be -1/(M*(M-2)) , otherwise a contradiction with special cases.

@Bendesarts   Do all things similarly, only replace  eq  by  lhs(eq_liaison[1, 1])  and etc.

@Carl Love  Visually  f[`3x`,`4y`]  is identical to  f[`3x, 4y`]

@yihezhi 

f[`{6x}`]; f[`{3x,4y}`];