WeBWorK Main Forum

$M->order_LR

$M->order_LR

by Nandor Sieben -
Number of replies: 10
Could someone explain how ->order_LR works. I am trying to use it to find the rank but something is not working. Is there a better way to find the rank? The example code below produces 4 for both $r and for $rLR but the rank of the matrix should be 3 according to the Matlab clone Freemat:

  rank([0,3,2,0;-4,2,-10,7;2,0,6,-4;2,-3,2,-1])
  ans = 3 

I am testing this on version 2.7 .

DOCUMENT();
loadMacros(
"PGstandard.pl",
"MathObjects.pl",
"AnswerFormatHelp.pl",
"PGcourse.pl",
);
$showPartialCorrectAnswers = 1;
TEXT(beginproblem());

Context('Matrix');

$aa=Matrix([[0,3,2,0],[-4,2,-10,7],[2,0,6,-4],[2,-3,2,-1]]);
$r=$aa->order_LR;

$aaLR=$aa->decompose_LR;
$rLR=$aaLR->order_LR;

Context()->texStrings;
BEGIN_TEXT
\( $aa, $r, $rLR \)

END_TEXT
Context()->normalStrings;

ENDDOCUMENT();
In reply to Nandor Sieben

Re: $M->order_LR

by Davide Cervone -
The problem is with pg/lib/MatrixReal1.pm, in its order_LR method (the MathObject Matrix simply calls the one from MatrixReal1). That routine looks at the diagonal of the R matrix from the LR decomposition and checks (from bottom up) if they are zero. But that check is explicitly against 0, rather than a fuzzy check for near zero. In the case of the matrix you have given, the (4,4) entry in R is 6.93889390390723E-17, which is not exactly zero, so order_LR counts that as non-zero and gives you order 4 rather than order 3.

I think that check should be a fuzzy check against the $zeroLevelTol. Something like 

    last ZERO if (CORE::abs($matrix->[0][$order][$order]) >= $main::zeroLevelTol);
Davide
In reply to Davide Cervone

Re: $M->order_LR

by Nandor Sieben -
What is the workaround? Would a call to $M->normalize(nullvector) help? Perhaps using rref($A) of MatrixReduce.pl with fraction context? Could we add a $M->rank function to Matrix math objects that works reliably with matrices containing small integer entries.
In reply to Nandor Sieben

Re: $M->order_LR, rank

by Nandor Sieben -
I have an attempt for a rank function. Davide, could you add something like this to the matrix mathobjects.

DOCUMENT();
loadMacros(
"PGstandard.pl",
"MathObjects.pl",
"AnswerFormatHelp.pl",
"PGcourse.pl",
"MatrixReduce.pl",
);
$showPartialCorrectAnswers = 1;
TEXT(beginproblem());

Context('Matrix');

sub rank {
my $M = shift;
Context('Fraction');
$M = apply_fraction_to_matrix_entries($M);
my $Mred=rref($M);
Context('Matrix');
my ($row,$col) = $Mred->dimensions;
{
  my $sum=0;
  for (my $i=1; $i<=$col; $i++) {
    $sum += abs($Mred->element($row,$i));
  }
  if (!$sum) {
    $row--;
    redo;
  }
}
$row;
}

$aa=Matrix([[0,3,2,0],[-4,2,-10,7],[2,0,6,-4],[2,-3,2,-1]]);
$r=$aa->order_LR;
$ra = rank($aa);

Context()->texStrings;
BEGIN_TEXT
\( $aa,$r,$ra \)
END_TEXT

Context()->normalStrings;

ENDDOCUMENT();
In reply to Nandor Sieben

Re: $M->order_LR, rank

by Nandor Sieben -
Chris Heckman sent me this alternate implementation.

sub rank {
   my $R = rref(shift);
   my ($m, $n) = $R-> dimensions;
   my $r = 1; 
   my $c = 1;
   while (($r <= $m) && ($c <= $n)) {
      $r++ if ($R->element($r,$c) > 0.5);
      $c ++;
      }
   $r - 1;
}

In reply to Nandor Sieben

Re: $M->order_LR

by Davide Cervone -
What is the workaround?

I think the work-around is to edit pg/lib/MatrixReal1.pm as I indicated above. It would be possible to write alternative routines, as you have suggested, but I think it is better to fix the existing one.
In reply to Davide Cervone

Re: $M->order_LR

by Nathan Wallach -
I wanted to test the rank of a matrix without the need to assume the values can be converted to fractions, as necessary in the code proposed by Nandor Sieben.

I also ran into a limitation that order_LR apparently (even with the fuzzy comparison to 0) only works well for square matrices.

The following is a provisional version of rank.pl which seem to work for me in my first test cases.

It has header code to handle the case of non-square matrices by adding 0 rows when there are more rows than columns. (The case of more columns than rows is handled by calling rank() recursively on the transpose matrix).

The indexing in the code is done using MathObjects 1-based indexing for rows/columns, so somewhat different that what was in MatrixReal1.pm which used Perl's 0-based indexing.

Testing and comments on this code would be appreciated.

From prior experience with Matlab, which allows passing in a tolerance level as a second argument to rank(A,tol) it may be necessary to consider adding such an optional second argument which could adjust the values of zeroLevelTol for the Context in which the comparison to Real("0") is made.




# Coded by Nathan Wallach, May 2019
# based on Davide Cervone's recommendation from
# http://webwork.maa.org/moodle/mod/forum/discuss.php?d=3194
# as the current $Matrix->order_LR does not use fuzzy comparisons
# so does not give good results.

sub rank { # It assumes that a MathObject matrix object is sent
my $MM1 = shift;

if ( $MM1->class ne "Matrix") {
return -1; # Not a matrix
}

# For the case it is square
my $MM = $MM1;

my ($Rrows,$Rcols) = $MM1->dimensions;

if ( ( $Rrows <= 0 ) || ( $Rcols <= 0 ) ) {
return -1; # Not a matrix
}

if ( $Rrows < $Rcols ) {
# pad to make it square
my @rows = ();
my $i = 1;
for ( $i = 1 ; $i <= $Rrows ; $i++ ) {
push( @rows, $MM1->row($i) );
}
while ( $i <= $Rrows ) {
# pad with zero rows
push( @rows, ( 0 * $MM1->row(1) ) );
$i++;
}
$MM = Matrix( @rows );
} elsif ( $Rrows > $Rcols ) {
return( rank( $MM1->transpose ) );
}

# Davide's approach from http://webwork.maa.org/moodle/mod/forum/discuss.php?d=3194
my $tempR = $MM->R;
($Rrows,$Rcols) = $tempR->dimensions;
my $rank;

for ( $rank = $Rrows ; $rank >= 1; $rank-- ) {
last if ( $tempR->element($rank,$rank) != Real("0") );
}
return( $rank );
}

1;
In reply to Nathan Wallach

Re: $M->order_LR

by Nathan Wallach -
See: https://github.com/openwebwork/pg/issues/409

Additional testing / feedback for the proposed rank and basis grading code would be appreciated.
In reply to Nathan Wallach

Re: $M->order_LR

by Ping-Shun Chan -

For the problem Library/NAU/setLinearAlgebra/expandBasis.pg

Under WeBWorK 2.15, and with random seed 1253, I get the error "Your answer isn't a matrix, it looks like a list of lists of numbers".

Somehow replacing the default rank.pl with Nathan's rank.pl fixes the problem.

(Thank you, Nathan!)