Double transpose on right_kernel()

[Rob Beezer]

I have been walking through some of the matrix code before adding some
enhancements.

A right_kernel() will take the transpose of a matrix and call
left_kernel(). left_kernel() (or kernel()) generally takes a
transpose before calling some other code to row-reduce the matrix, and
so on. It appears that most of the specialized code actually
constructs a right kernel: IML for dense integer/rational matrices,
PARI for those with entries in number fields and original Sage code
for other fields. The effect is that a matrix is often (always?)
transposed twice if you are really after a right kernel.

I built a 2500 x 2500 matrix and multiplied it by its transpose so I
could work on a symmetric matrix (thus a left or right kernel will be
identical computations).

sage: a=random_matrix(ZZ,2500,x=10)
sage: m=a*a.transpose().change_ring(QQ)

sage: time m.left_kernel() --> 5.84 seconds
sage: time m.right_kernel() --> 8.58 seconds

Timings on just a simple transpose of a matrix of this size pretty
much account for all of the difference. It should therefore be
possible to get a right kernel in this case at around 3 seconds if two
transposes were eliminated. Furthermore, I've found all the
transposes a bit confusing as I chase my way up and down the class
hierarchy.

Proposal: Place the guts of kernel computations for each
(specialized) class into a right_kernel() method, where it would seem
to naturally belong. Mostly these would replace kernel() methods that
are computing left kernels. A call to kernel() or left_kernel()
should arrive at the top of the hierarchy where it would take a
transpose and call the (specialized) right_kernel(). So there
wouldn't be a change in behavior in routines currently calling kernel
() or left_kernel(), and we would retain Sage's preference for the
left by having the vanilla kernel() give back a left kernel. The
changes would be more efficient computationally and also make the code
more expressive of what is really happening when and where.

Comments, suggestions, warnings?

[David Joyner]

On Sun, Feb 22, 2009 at 1:30 AM, Rob Beezer <goo...@beezer.cotse.net> wrote:
>
> I have been walking through some of the matrix code before adding some
> enhancements.
>

...

>
> Proposal: Place the guts of kernel computations for each
> (specialized) class into a right_kernel() method, where it would seem
> to naturally belong. Mostly these would replace kernel() methods that

+1

> are computing left kernels. A call to kernel() or left_kernel()
> should arrive at the top of the hierarchy where it would take a
> transpose and call the (specialized) right_kernel(). So there
> wouldn't be a change in behavior in routines currently calling kernel
> () or left_kernel(), and we would retain Sage's preference for the
> left by having the vanilla kernel() give back a left kernel. The
> changes would be more efficient computationally and also make the code
> more expressive of what is really happening when and where.
>
> Comments, suggestions, warnings?

(1) Add very detailed comments to the docstring describing exactly
what is going on.

(2) Test before and after behavior/timings to be sure the patch does as you
hope.


>
> >
>

[Rob Beezer]

The transposes have been sped-up nicely by Yann Laigle-Chapuy
http://trac.sagemath.org/sage_trac/ticket/5345
http://trac.sagemath.org/sage_trac/ticket/5369

I've added the refactoring into Trac as #5381
http://trac.sagemath.org/sage_trac/ticket/5381

Rob

[David Kohel]

Hi,

As I have been creditted by William with the free modules
as row vectors and default left kernel for matrices (under
the influence of Magma), let me propose a coherent way
to relax this design decision.

It should be easy to add an argument to free modules to
create them as column spaces and give a nice latex and
(not so nice) ascii art printing. Thus right_kernel should
return such a module or vector space.

Arbitrary matrices will then have both left and right kernels,
and the composition will be well-defined. There has been
a huge amount of effort going into calculus to make Sage
suitable (and a natural choice) for undergraduate calculus
teaching.

A much smaller amount of effort would make it suitable to
map any undergraduate linear algebra textbook exercise
in terms of row OR column vectors to a Sage exercise
without convoluted applications of transposes.

Except for details like passing the column/row orientation
to subspaces and caching issues, the main changes are
largely restricted to __repr__ and __mul__ (for compatibility
with matrix multiplication). If someone has two days to
spare (one for implementation and one for documentation
and bug tracking), then I would highly recommend this
complementary project to the right_kernel cleanup.

There would be numerous advantages also to research,
e.g. in the theory of theta functions or Siegel modular
forms it would be a huge pain to correctly transpose the
standard column vector notation into reliable code.

--David

[Rob Beezer]

Hi David,

Thanks for your comments. I'm very interested in making Sage more
approachable for students, so I'm all for anything that improves that,
without hampering research applications.

In my linear algebra textbook, I've been dogmatic about every vector
being a column vector, even going so far as constructing row spaces as
sets of column vectors. Definitely over the top, and arguably wrong.
But I've been trying to spare the students some initial confusion
about rows and columns and when to transpose, and when not to.

Would it be going too far to suggest vectors could be tagged as "row"
or "column" (or None)? A right kernel would return a module/subspace
with basis vectors that were tagged as columns and a left kernel would
return a module/subspace with basis vectors that were tagged as rows.
A column space would have columns as basis vectors and a row space
would have basis vectors that are rows. I could then see sanity
checks other places, like the product of a matrix with a vector on the
right or left could check for the vector being a column or a row. An
inner product could require a row vector, then a column vector, in
that order as arguments. A transpose method on vectors would just
toggle the tag. And so on. These checks could be phased in, choosing
where to query the tag and where to respect the None tag. I could see
a payoff in preventing some incompatible operations in user-level code
with clear messages (eg "attempting to multiply a column vector on the
left of a matrix").

Would this be too disruptive in other places?

Printing. In my cynical moments, I've long thought row vectors are
just a typographical convenience, and certainly the (1, 2, 3)^t
construction for a printed column vector must be a publisher's
preference. But here it might make sense. Just print vectors tagged
"column" exactly as is done now, but with "^t" tacked on?

I hope I've understood your comments correctly and haven't run off too
far in the wrong direction.

Rob

[William Stein]

It was my impression that your proposal above is identical to the one
that David Kohel just proposed.

> I hope I've understood your comments correctly and haven't run off too
> far in the wrong direction.

I think you have.