nullspace computation
Shunsuke Tsuchioka
I am a graduate student in the math department of Kyoto university,
and want to calculate a basis of the nullspace
\{ x\in F^n \mid Ax=0 \}
for a given m times n matrix A whose entries are in a field F.
What is a good way to calculate a nullspace basis
for a given matrix A by LinBox?
The following is what I have tried:
Fortunately in my research, a nullspace is mathematically guaranteed
to be 1-dimensional. Thus, I first thought "linbox/solutions/solve.h"
may help. But the code like the following
--------------------------------------------------------------
typedef LinBox::Modular< double > Field;
Field F(q);
LinBox::DenseMatrix< Field > A(F,row,col);
std::vector<Field::Element> x( A.coldim() ),b( A.rowdim() );
A.setEntry( abracadabra, xyz, ... )
for( std::vector<Field::Element>::iterator it=b.begin(); it != b.end(); ++it ) { *it = 0; }
LinBox::solve( x, A, b, LinBox::Method::BlasElimination() );
---------------------------------------------------------------
always returns the solution x = 0 (of course, it is a valid solution).
I am sorry for such a basic question, but
it would be great if you colud help me.
Shunsuke, TSUCHIOKA
Research Institute for Mathematical Sciences Kyoto University
William Stein
>
> Dear Linbox Support team,
>
> I am a graduate student in the math department of Kyoto university,
> and want to calculate a basis of the nullspace
>
> \{ x\in F^n \mid Ax=0 \}
>
> for a given m times n matrix A whose entries are in a field F.
> What is a good way to calculate a nullspace basis
> for a given matrix A by LinBox?
>
You might want to try doing this with Sage (http://sagemath.org), which
is a free open source program with a very nice interpreter interface that
is built on top of Linbox and provides the above functionality. E.g.,
dhcp46-76:hnf was$ sage
----------------------------------------------------------------------
| SAGE Version 2.10.1, Release Date: 2008-02-02 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
sage: a = random_matrix(QQ, 5,3); a
[ -1 -2 1]
[-1/2 0 0]
[ 1/2 -2 1]
[ 1 2 1]
[ 0 -2 -1/2]
sage: a.kernel()
Vector space of degree 5 and dimension 2 over Rational Field
Basis matrix:
[ 1 0 -8/5 9/5 12/5]
[ 0 1 -1/5 3/5 4/5]
sage: a.kernel().basis()
[
(1, 0, -8/5, 9/5, 12/5),
(0, 1, -1/5, 3/5, 4/5)
]
Note that kernel returns the set of row vectors v such that v*A = 0.
You can do quite large problems using Sage (which again, builds on
Linbox, so it's
OK that I'm posting here):
sage: a = random_matrix(QQ, 100, 102)
sage: time v = a.kernel()
CPU times: user 0.04 s, sys: 0.02 s, total: 0.06 s
sage: a = random_matrix(QQ, 400, 402)
sage: time v = a.kernel()
CPU times: user 0.59 s, sys: 0.09 s, total: 0.68 s
Working modulo 97:
sage: a = random_matrix(GF(97), 100, 102)
sage: time v = a.kernel()
CPU times: user 0.04 s, sys: 0.00 s, total: 0.04 s
Wall time: 0.06
If the specific functionality (the field etc) that you need isn't available
in Sage or is too slow, let me know, and Clement Pernet (one of the main Linbox
developers) and I can very likely add support for using that
functionality of Linbox in Sage.
> The following is what I have tried:
>
> Fortunately in my research, a nullspace is mathematically guaranteed
> to be 1-dimensional. Thus, I first thought "linbox/solutions/solve.h"
> may help. But the code like the following
>
> --------------------------------------------------------------
> typedef LinBox::Modular< double > Field;
> Field F(q);
> LinBox::DenseMatrix< Field > A(F,row,col);
> std::vector<Field::Element> x( A.coldim() ),b( A.rowdim() );
> A.setEntry( abracadabra, xyz, ... )
> for( std::vector<Field::Element>::iterator it=b.begin(); it != b.end(); ++it ) { *it = 0; }
> LinBox::solve( x, A, b, LinBox::Method::BlasElimination() );
> ---------------------------------------------------------------
>
> always returns the solution x = 0 (of course, it is a valid solution).
>
> I am sorry for such a basic question, but
> it would be great if you colud help me.
>
> Shunsuke, TSUCHIOKA
> Research Institute for Mathematical Sciences Kyoto University
>
> >
>
--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
Shunsuke Tsuchioka
Thank you for your help.
SAGE project sounds good as you wrote to me and I will try.
But, if possible, is there a way to calculate nullspace basis by LinBox C++ code?
Because my research program requires many headers (*.hpp) which I wrote about a year ago,
it seems to good to continue writing in C++ language.
Thank you again.
B.David Saunders
It is something we should provide directly and I think will soon do so.
In the meantime, for your problem you could strip off the last column, A
= ( A' b ), and solve A' x = b. Then x with -1 adjoined is a nullspace
vector. That will work as long as column b is in the span of the other
columns. Failing that you could solve Ax = Ay, where y is a random
vector. For your nullity = 1 matrix, it is not likely that linbox will
return x being equal to y. Thus with high probability x-y will be a
nontrivial nullspace vector.
I mention that SAGE is a wonderful tool. You can call C++ code from
SAGE and SAGE code from C++ so it is not incompatible with using your
existing code.
-dave
David Saunders
Dept. of Computer and Information Sciences, Univ. of Delaware
email: saun...@cis.udel.edu Phone: 302-831-6238
B.David Saunders
For the dense case it is a no-brainer that we should provide nullspace
basis (even though its entries may be huge). For the sparse case the
basis can be very much larger than the input and some user problems can
be solved by SAMPLING the nullspace rather than first computing a
complete basis. As that can be much more efficient of time and memory,
we should be sure we provide for that (and encourage it). Here is my
proposal for functions in solutions/ for solving and for nullspace:
1. Special function for nonsingular matrices.
solveNonsingular(x, A, b) // compute x = A^{-1}b.
2. General case, where n is the (right) nullity of matrix/blackbox A.
solve(x, A, b) // compute ARBITRARY x, such that Ax = b. If no such
x exists, set x = 0.
Here 'arbitrary' means easiest for the implementor. This is what we
offer now and what returns x = 0 when b = 0 and the method is
blaselimination.
solve(x, A, b, M) // user can use the method M to specify other
handling of the inconsistent case (eg to give moore-penrose solution or
to indicate inconsistency in the M on output - ignoring x).
nullspaceRandom(x, A) // compute a RANDOM element of nullspace(A).
Whatever else 'random' means, it implies that, with high likelihood, n
calls will produce a nullspace basis.
nullspaceBasis(X, A) // compute a n column, full column rank
DenseMatrix X such that AX = 0. (with BlasBlackbox becoming the
DenseMatrix soon, and, for the case n = 0, with an X having zero columns
being allowed by the representation!) Note that this implies a memory
model in which the function is allowed to resize the output X. ...and
we also give suitable remarks to the user in the comments about the case
of a sparse matrix with large nullity ... or the case of a dense integer
matrix with large nullity where again the basis can require n times more
storage than A.
Note that solveNonsingular allows for some extra efficiency in case
nonsingularity is known a priori. The user's assertion of
nonsingularity could be embedded in the method. But I think it is good
to embed it in the function name for clarity, ease of use with default
method, etc.
The three general functions solve(), nullspaceVector(), nullspaceBasis()
allow for good deployment of LinBox' various elimination based and
blackbox methods. (all these functions take the optional method
parameter of course, in the solutions/ directory style.) For the sparse
case it may be good - for efficiency, including block method efficiency
- to offer a version of nullspaceVector() that returns a clump of
(independent) nullspace vectors.
The following combination functions can then be left to the user or
offered as well.
solve(X, A, b) // compute X, a DenseMatrix of coldim n, such that Ax
= b, for each col x of X.
{ return X = (x + y_1, ..., x + y_n),
where the y_i are cols of nullspaceBasis(Y, A), and x is from
solve(x, A, b) // arbitrary
}
solveRandom(x, A, b) // compute a random element x of the solution
space (an affine plane). Remark: in some cases we have a better way
than combining an arbitrary solution and a random nullspace vector.
Note that this proposal relegates some other choices
-dave
PS. William, you are always welcome to discuss SAGE linear algebra
capability on this forum, without reservation or hesitation, whether
LinBox is involved or not. Out of curiosity, does your current
nullspace basis code use linbox at all (eg start with a linbox LU)?
PPS. I suppose this ought to be a linbox-devel post. But I think in
practice the distinction between linbox-use and linbox-devel is not
happening nor worth the effort.
Michael Abshoff
>
> LinBoxers,
<SNIP>
> PS. William, you are always welcome to discuss SAGE linear algebra
> capability on this forum, without reservation or hesitation, whether
> LinBox is involved or not. Out of curiosity, does your current
> nullspace basis code use linbox at all (eg start with a linbox LU)?
Hi Dave,
Sage currently uses IML's p-adics nullspace per default (also for ZZ), but
has fallback methods. The excellent performace of IML for nullspaces come
up this week during Sage Days 7, but Clement can probably fill in more of
the details about that since I worked on other issues.
Cheers,
Michael
William Stein
>
> LinBoxers,
>
>
> PS. William, you are always welcome to discuss SAGE linear algebra
> capability on this forum, without reservation or hesitation, whether
> LinBox is involved or not.
Thanks! I'll do so :-)
> Out of curiosity, does your current
> nullspace basis code use linbox at all (eg start with a linbox LU)?
It depends on the base ring.
ZZ -- currently it uses PARI, since PARI has the best Hermite Normal
Form algorithm right now (of course that's changing soon!)
QQ -- clear denominators and use the nullspaceMP command in IML.
This is the fastest implementation I know of for dense matrices,
especially with large coefficient sizes.
All other rings (e.g., GF(p))-- compute the reduced row echelon form and read
off the kernel. Echelon form, e.g., for GF(p) is computed using Linbox.
In other cases they are computed using a block strassen algorithm that David
Harvard and Robert Bradshaw designed and implemented.
Thus for the specific question that started this thread, i.e., computing kernels
over GF(p), Sage *is* using Linbox to compute the echelon form, then at the end
reading off the kernel from that.
> PPS. I suppose this ought to be a linbox-devel post. But I think in
> practice the distinction between linbox-use and linbox-devel is not
> happening nor worth the effort.
If the traffic picks up a lot it would make sense to worry about it...
Clement Pernet
Please answer to linbox-devel only.
> 1. Special function for nonsingular matrices.
> solveNonsingular(x, A, b) // compute x = A^{-1}b.
>
Ok why not, although I am not sure there is much efficiency to gain from
the knowledge of the non-singularity.
> 2. General case, where n is the (right) nullity of matrix/blackbox A.
>
> solve(x, A, b) // compute ARBITRARY x, such that Ax = b. If no such
> x exists, set x = 0.
>
I disagree. The value of x should not be used as an exception flag.
Let's use exceptions or dedicated flags.
I would prefer exceptions: we have them, and the cost of catching them
is really not an issue.
> Here 'arbitrary' means easiest for the implementor. This is what we
> offer now and what returns x = 0 when b = 0 and the method is
> blaselimination.
>
> solve(x, A, b, M) // user can use the method M to specify other
> handling of the inconsistent case (eg to give moore-penrose solution or
> to indicate inconsistency in the M on output - ignoring x).
>
>
Ok with that.
> nullspaceRandom(x, A) // compute a RANDOM element of nullspace(A).
> Whatever else 'random' means, it implies that, with high likelihood, n
> calls will produce a nullspace basis.
>
> nullspaceBasis(X, A) // compute a n column, full column rank
> DenseMatrix X such that AX = 0. (with BlasBlackbox becoming the
> DenseMatrix soon, and, for the case n = 0, with an X having zero columns
> being allowed by the representation!) Note that this implies a memory
> model in which the function is allowed to resize the output X. ...and
> we also give suitable remarks to the user in the comments about the case
> of a sparse matrix with large nullity ... or the case of a dense integer
> matrix with large nullity where again the basis can require n times more
> storage than A.
>
> Note that solveNonsingular allows for some extra efficiency in case
> nonsingularity is known a priori. The user's assertion of
> nonsingularity could be embedded in the method. But I think it is good
> to embed it in the function name for clarity, ease of use with default
> method, etc.
>
Non-singularity only implies simpler code ; but I doubt it enables
better efficiency. Correct me if I am wrong.
> The three general functions solve(), nullspaceVector(), nullspaceBasis()
> allow for good deployment of LinBox' various elimination based and
> blackbox methods. (all these functions take the optional method
> parameter of course, in the solutions/ directory style.) For the sparse
> case it may be good - for efficiency, including block method efficiency
> - to offer a version of nullspaceVector() that returns a clump of
> (independent) nullspace vectors.
>
> The following combination functions can then be left to the user or
> offered as well.
>
> solve(X, A, b) // compute X, a DenseMatrix of coldim n, such that Ax
> = b, for each col x of X.
> { return X = (x + y_1, ..., x + y_n),
> where the y_i are cols of nullspaceBasis(Y, A), and x is from
> solve(x, A, b) // arbitrary
> }
>
>
Yes, let's provide it.
btw, we also have to provide the solve (X,A,B) where X and B are n*k
matrices (solving for multiple columns).
> solveRandom(x, A, b) // compute a random element x of the solution
> space (an affine plane). Remark: in some cases we have a better way
> than combining an arbitrary solution and a random nullspace vector.
>
> Note that this proposal relegates some other choices
>
My proposal in the next post (to avoid spaghetti message).
Clément