spral_scaling - Sparse matrix scalings

This package generates various scalings (and matchings) of real sparse matrices.

Given a symmetric matrix \(A\), it finds a diagonal matrix \(D\) such that the scaled matrix

\[\hat{A} = DAD\]

has specific numerical properties.

Given an unsymmetric or rectangular matrix \(A\), it finds diagonal matrices \(D_r\) and \(D_c\) such that the scaled matrix

\[\hat{A} = D_r A D_c\]

has specific numerical properties.

The specific numerical properties delivered depends on the algorithm used:

Matching-based

algorithms scale \(A\) such that the maximum (absolute) value in each row and column of \(\hat{A}\) is exactly \(1.0\), where the entries of maximum value form a maximum cardinality matching. The Hungarian algorithm delivers an optimal matching slowly, whereas the auction algorithm delivers an approximate matching quickly.

Norm-equilibration

algorithms scale \(A\) such that the infinity norm of each row and column of \(\hat{A}\) is \(1.0\pm \tau\) (for some user specified tolerance \(\tau\)).

Auction Algorithm

Routines

subroutine  auction_scale_sym(n, ptr, row, val, scaling, options, inform[, match])

Find a matching-based symmetric scaling using the auction algorithm.

The scaled matrix is such that the entry of maximum absolute value in each row and column is (approximately) \(1.0\).

Parameters:

• n [integer ,in] :: number of columns in \(A\).

• ptr (n+1) [integer ,in] :: columns pointers for \(A\) (see CSC format).

• row (ptr(n+1)-1) [integer ,in] :: row indices for \(A\) (see CSC format).

• val (ptr(n+1)-1) [real ,in] :: non-zero values for \(A\) (see CSC format).

• scaling (n) [real ,out] :: returns scaling found by routine.

• options [auction_options ,in] :: controls behaviour of routine.

• inform [auction_inform ,out] :: returns information on execution of routine.

Options:

match (n) [integer ,out] :: returns matching found by routine. Row i is matched to column match(i), or is unmatched if match(i)=0.

subroutine  auction_scale_unsym(m, n, ptr, row, val, rscaling, cscaling, options, inform[, match])

Find a matching-based unsymmetric scaling using the auction algorithm.

The scaled matrix is such that the entry of maximum absolute value in each row and column is (approximately) \(1.0\).

Parameters:

• m [integer ,in] :: number of rows in \(A\)

• n [integer ,in] :: number of columns in \(A\)

• ptr (n+1) [integer ,in] :: columns pointers for \(A\) (see CSC format)

• row (ptr(n+1)-1) [integer ,in] :: row indices for \(A\) (see CSC format)

• val (ptr(n+1)-1) [real ,in] :: non-zero values for \(A\) (see CSC format)

• rscaling (m) [real ,out] :: returns row scaling found by routine

• cscaling (n) [real ,out] :: returns column scaling found by routine

• options [auction_options ,in] :: controls behaviour of routine

• inform [auction_inform ,out] :: returns information on execution of routine

Options:

match (n) [integer ,out] :: returns matching found by routine. Row i is matched to column match(i), or is unmatched if match(i)=0.

Data-types

type  auction_options

Used to specify options to the routines auction_scale_sym() and auction_scale_unsym(). Components are automatically given default values in the definition of the type.

Please refer to the method section for details on how these parameters are used.

Type fields:

• % eps_initial [real ,default=0.01] :: initial value of improvement parameter \(\epsilon\).

• % max_iterations [integer ,default=30000] :: maximum number of iterations.

• % max_unchanged (3) [integer ,default=10,100,100] :: together with min_proportion(:), specifies termination conditions.

• % min_proportion (3) [real ,default=0.9,0.0,0.0] :: together with max_unchanged(:), specifies termination conditions.

type  auction_inform

Used to return information about the execution of the algorithm.

Type fields:

• % flag [integer ] :: gives the exit status of the algorithm (see table below)

• % iterations [integer ] :: number of iteration performed.

• % matched [integer ] :: number of rows and columns that have been matched.

• % stat [integer ] :: holds the Fortran stat parameter in the event of an allocation failure (set to 0 otherwise).

• % unmatchable :: holds the number of columns designated as unmatchable (there is no way to match it that improves the quality of the matching).

Note: As the algorithm may terminate before a full matching is obtained, inform%matched provides only a lower bound on the structural rank. However, inform%unmatchable provides an approximate lower bound on the structural rank deficiency.

inform%flag	Return status
0	Success.
-1	Allocation error. Fortran stat value is returned in inform%stat.

Example

The following code shows an example usage of auction_scale_sym().

! examples/Fortran/scaling/auction_sym.f90 - Example code for SPRAL_SCALING
program auction_scale_sym_example
   use spral_scaling
   use spral_matrix_util, only : print_matrix, &
                                 SPRAL_MATRIX_REAL_SYM_INDEF
   implicit none

   ! Derived types
   type (auction_options)  :: options
   type (auction_inform)   :: inform

   ! Parameters
   integer, parameter :: wp = kind(0.0d0)

   ! Matrix data
   integer :: n, ptr(6), row(8)
   real(wp) :: val(8)

   ! Other variables
   integer :: match(5), i, j
   real(wp) :: scaling(5)

   ! Data for symmetric matrix:
   ! ( 2  1         )
   ! ( 1  4  1    8 )
   ! (    1  3  2   )
   ! (       2      )
   ! (    8       2 )
   n = 5
   ptr(1:n+1)        = (/ 1,        3,             6,      8,8,   9 /)
   row(1:ptr(n+1)-1) = (/ 1,   2,   2,   3,   5,   3,   4,   5   /)
   val(1:ptr(n+1)-1) = (/ 2.0, 1.0, 4.0, 1.0, 8.0, 3.0, 2.0, 2.0 /)
   write(*, "(a)") "Initial matrix:"
   call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)

   ! Perform symmetric scaling
   call auction_scale_sym(n, ptr, row, val, scaling, options, inform, &
      match=match)
   if(inform%flag<0) then
      write(*, "(a, i5)") "auction_scale_sym() returned with error ", &
         inform%flag
      stop
   endif

   ! Print scaling and matching
   write(*,"(a,10i10)")    'Matching:', match(1:n)
   write(*,"(a,10es10.2)") 'Scaling: ', scaling(1:n)

   ! Calculate scaled matrix and print it
   do i = 1, n
      do j = ptr(i), ptr(i+1)-1
         val(j) = scaling(i) * val(j) * scaling(row(j))
      end do
   end do
   write(*, "(a)") "Scaled matrix:"
   call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)

end program auction_scale_sym_example

The above code produces the following output:

Initial matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
1:   2.0000E+00   1.0000E+00
2:   1.0000E+00   4.0000E+00   1.0000E+00                8.0000E+00
3:                1.0000E+00   3.0000E+00   2.0000E+00
4:                             2.0000E+00
5:                8.0000E+00                             2.0000E+00
Matching:         1         5         4         3         2
Scaling:   7.07E-01  1.62E-01  2.78E-01  1.80E+00  7.72E-01
Scaled matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
1:   1.0000E+00   1.1443E-01
2:   1.1443E-01   1.0476E-01   4.5008E-02                1.0000E+00
3:                4.5008E-02   2.3204E-01   1.0000E+00
4:                             1.0000E+00
5:                1.0000E+00                             1.1932E+00

Method

This algorithm finds a fast approximation to the matching and scaling produced by the HSL package MC64. If an optimal matching is required, use the Hungarian algorithm instead. The algorithm works by solving the following maximum product optimization problem using an auction algorithm. The scaling is derived from the dual variables associated with the solution.

\[\begin{split}\max_{\sigma} & \prod_{i=1}^m\prod_{j=1}^n |a_{ij}|\sigma_{ij} & \\ \mathrm{s.t.} & \sum_{i=1}^m\sigma_{ij} = 1, & \forall j=1,n \\ & \sum_{j=1}^n\sigma_{ij} = 1, & \forall i=1,m \\ & \sigma_{ij} \in \{0,1\}.\end{split}\]

The array \(\sigma\) gives a matching of rows to columns.

By using the transformation

\[w_{ij} = \log c_j - \log |a_{ij}|,\]

where \(c_j = \max_i |a_{ij}|\), the maximum product problem in \(a_{ij}\) is replaced by a minimum sum problem in \(w_{ij}\) where all entries are positive. By standard optimization theory, there exist dual variables \(u\) and \(v\) corresponding to the constraints that satisfy the first order optimality conditions

\[\begin{split}w_{ij} - u_i - v_j = 0, & \mbox{ if } \sigma_{ij }=1, \\ w_{ij} - u_i - v_j \ge 0, & \mbox{ if } \sigma_{ij }=0.\end{split}\]

To obtain a scaling we define scaling matrices \(D_r\) and \(D_c\) as

\[ \begin{align}\begin{aligned}d^r_i = e^{u_i},\\d^c_i = e^{v_i}.\end{aligned}\end{align} \]

If a symmetric scaling is required, we average these as

\[d_i = \sqrt{d^r_id^c_i}.\]

By the first order optimality conditions, these scaling matrices guarantee that

\[\begin{split}d^r_i|a_{ij}|d^c_j = 1, && \mbox{if } \sigma_{ij}=1, \\ d^r_i|a_{ij}|d^c_j \le 1, && \mbox{if } \sigma_{ij}=0.\end{split}\]

To solve the minimum sum problem an auction algorithm is used. The algorithm is not guaranteed to find an optimal matching. However it can find an approximate matching very quickly. A matching is maintained along with the row pricing vector \(u\). In each major iteration, we loop over each column in turn. If the column \(j\) is unmatched, we calculate the value \(p_i = w_{ij} - u_i\) for each entry and find the maximum across the column. If this maximum is positive, the current matching can be improved by matching column \(j\) with row \(i\). This may mean that the previous match of row \(i\) now becomes unmatched. We update the price of row \(i\), that is \(u_i\), to reflect this new benefit and continue to the next column.

To prevent incremental shuffling, we insist that the value of adding a new column is at least a threshold value \(\epsilon\) above zero, where \(\epsilon\) is based on the last iteration in which row \(i\) changed its match. This is done by adding \(\epsilon\) to the price \(u_i\), where \(\epsilon = \mathrm{options\%eps\_initial} + \mathrm{itr} / (n+1)\), where itr is the current iteration number.

The algorithm terminates if any of the following are satsified:

• All entries are matched.

• The number of major iterations exceeds options%max_iterations.

• At least options%max_unchanged(1) iterations have passed without the cardinality of the matching increasing, and the proportion of matched columns is options%min_proportion(1).

• At least options%max_unchanged(2) iterations have passed without the cardinality of the matching increasing, and the proportion of matched columns is options%min_proportion(2).

• At least options%max_unchanged(3) iterations have passed without the cardinality of the matching increasing, and the proportion of matched columns is options%min_proportion(3).

The different combinations given by options%max_unchanged(1:3) and options%min_proportion(1:3) allow a wide range of termination heuristics to be specified by the user depending on their particular needs. Note that the matching and scaling produced will always be approximate as \(\epsilon\) is non-zero.

Further details are given in the following paper:

[1]

J.D. Hogg and J.A. Scott. (2014). On the efficient scaling of sparse symmetric matrices using an auction algorithm. RAL Technical Report RAL-P-2014-002. [STFC TR]

Norm-equilibration Algorithm

Routines

subroutine  equilib_scale_sym(n, ptr, row, val, scaling, options, inform)

Find a matching-based symmetric scaling using the norm-equilibration algorithm.

The scaled matrix is such that the infinity norm of each row and column are equal to \(1.0\).

Parameters:

• n [integer ,in] :: number of columns in \(A\).

• ptr (n+1) [integer ,in] :: columns pointers for \(A\) (see CSC format).

• row (ptr(n+1)-1) [integer ,in] :: row indices for \(A\) (see CSC format).

• val (ptr(n+1)-1) [real ,in] :: non-zero values for \(A\) (see CSC format).

• scaling (n) [real ,out] :: returns scaling found by routine.

• options [equilib_options ,in] :: controls behaviour of routine.

• inform [equilib_inform ,out] :: returns information on execution of routine.

subroutine  equilib_scale_unsym(m, n, ptr, row, val, rscaling, cscaling, options, inform)

Find a matching-based unsymmetric scaling using the norm-equilibration algorithm.

The scaled matrix is such that the infinity norm of each row and column are equal to \(1.0\).

Parameters:

• m [integer ,in] :: number of rows in \(A\).

• n [integer ,in] :: number of columns in \(A\).

• ptr (n+1) [integer ,in] :: columns pointers for \(A\) (see CSC format).

• row (ptr(n+1)-1) [integer ,in] :: row indices for \(A\) (see CSC format).

• val (ptr(n+1)-1) [real ,in] :: non-zero values for \(A\) (see CSC format).

• rscaling (m) [real ,out] :: returns row scaling found by routine.

• cscaling (n) [real ,out] :: returns column scaling found by routine.

• options [equilib_options ,in] :: controls behaviour of routine.

• inform [equilib_inform ,out] :: returns information on execution of routine.

Data-types

type  equilib_options

Used to specify options to the routines equilib_scale_sym() and equilib_scale_unsym(). Components are automatically given default values in the definition of the type.

Please refer to the method section for details on how these parameters are used.

Type fields:

• % max_iterations [integer ,default=10] :: maximum number of iterations.

• % tol [real ,default=1e-8] :: convergence tolerance \(tau\).

type  equilib_inform

Used to return information about the execution of the algorithm.

Type fields:

• % flag [integer ] :: gives the exit status of the algorithm (see table below)

• % iterations [integer ] :: number of iteration performed.

• % stat [integer ] :: holds the Fortran stat parameter in the event of an allocation failure (set to 0 otherwise).

inform%flag	Return status
0	Success.
-1	Allocation error. Fortran stat value is returned in inform%stat.

Example

The following code shows an example usage of equilib_scale_sym().

! examples/Fortran/scaling/equilib_sym.f90 - Example code for SPRAL_SCALING
program equilib_scale_sym_example
   use spral_scaling
   use spral_matrix_util, only : print_matrix, &
                                 SPRAL_MATRIX_REAL_SYM_INDEF
   implicit none

   ! Derived types
   type (equilib_options)  :: options
   type (equilib_inform)   :: inform

   ! Parameters
   integer, parameter :: wp = kind(0.0d0)

   ! Matrix data
   integer :: n, ptr(6), row(8)
   real(wp) :: val(8)

   ! Other variables
   integer :: i, j
   real(wp) :: scaling(5)

   ! Data for symmetric matrix:
   ! ( 2  1         )
   ! ( 1  4  1    8 )
   ! (    1  3  2   )
   ! (       2      )
   ! (    8       2 )
   n = 5
   ptr(1:n+1)        = (/ 1,        3,             6,      8,8,   9 /)
   row(1:ptr(n+1)-1) = (/ 1,   2,   2,   3,   5,   3,   4,   5   /)
   val(1:ptr(n+1)-1) = (/ 2.0, 1.0, 4.0, 1.0, 8.0, 3.0, 2.0, 2.0 /)
   write(*, "(a)") "Initial matrix:"
   call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)

   ! Perform symmetric scaling
   call equilib_scale_sym(n, ptr, row, val, scaling, options, inform)
   if(inform%flag<0) then
      write(*, "(a, i5)") "equilib_scale_sym() returned with error ", &
         inform%flag
      stop
   endif

   ! Print scaling and matching
   write(*,"(a,10es10.2)") 'Scaling: ', scaling(1:n)

   ! Calculate scaled matrix and print it
   do i = 1, n
      do j = ptr(i), ptr(i+1)-1
         val(j) = scaling(i) * val(j) * scaling(row(j))
      end do
   end do
   write(*, "(a)") "Scaled matrix:"
   call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)

end program equilib_scale_sym_example

The above code produces the following output:

Initial matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
1:   2.0000E+00   1.0000E+00
2:   1.0000E+00   4.0000E+00   1.0000E+00                8.0000E+00
3:                1.0000E+00   3.0000E+00   2.0000E+00
4:                             2.0000E+00
5:                8.0000E+00                             2.0000E+00
Scaling:   7.07E-01  3.54E-01  5.77E-01  8.66E-01  3.54E-01
Scaled matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
1:   1.0000E+00   2.5000E-01
2:   2.5000E-01   5.0000E-01   2.0412E-01                1.0000E+00
3:                2.0412E-01   1.0000E+00   9.9960E-01
4:                             9.9960E-01
5:                1.0000E+00                             2.5000E-01

Method

This algorithm is very similar to that used by the HSL routine MC77. An iterative method is used to scale the infinity norm of both rows and columns to \(1.0\) with an asymptotic linear rate of convergence of \(\frac{1}{2}\), preserving symmetry if the matrix is symmetric.

For unsymmetric matrices, the algorithm outline is as follows:

• \(D_r^{(0)} = I, D_c^{(0)}=I\)

• for (\(k=1,\) options%max_iterations)

  • \(A^{(k-1)} = D_r^{(k-1)} A D_c^{(k-1)}\)

  • \((D_r^{(k)})_{ii} = (D_r^{(k-1)})_{ii}\; /\; \sqrt{\max_j(A^{(k-1)})_{ij}}\)

  • \((D_c^{(k)})_{jj} = (D_c^{(k-1)})_{jj}\; /\; \sqrt{\max_i(A^{(k-1)})_{ij}}\)

  • if ( \(|1-\|A^{(k-1)}\|_{\max}|\le\) options%tol ) exit

For symmetric matrices, \(A^{(k-1)}\) is symmetric, so \(D_r^{(k)} = D_c^{(k)}\), and some operations can be skipped.

Further details are given in the following paper:

[2]

16. Knight, D. Ruiz and B. Ucar. (2012). A symmetry preserving algorithm for matrix scaling. INRIA Research Report 7552. [INRIA TR]

Hungarian Algorithm

Routines

subroutine  hungarian_scale_sym(n, ptr, row, val, scaling, options, inform[, match])

Find a matching-based symmetric scaling using the Hungarian algorithm.

The scaled matrix is such that the entry of maximum absolute value in each row and column is \(1.0\).

Parameters:

• n [integer ,in] :: number of columns in \(A\).

• ptr (n+1) [integer ,in] :: columns pointers for \(A\) (see CSC format).

• row (ptr(n+1)-1) [integer ,in] :: row indices for \(A\) (see CSC format).

• val (ptr(n+1)-1) [real ,in] :: non-zero values for \(A\) (see CSC format).

• scaling (n) [real ,out] :: returns scaling found by routine.

• options [hungarian_options ,in] :: controls behaviour of routine.

• inform [hungarian_inform ,out] :: returns information on execution of routine.

Options:

match (n) [integer ,out] :: returns matching found by routine. Row i is matched to column match(i), or is unmatched if match(i)=0.

subroutine  hungarian_scale_unsym(m, n, ptr, row, val, rscaling, cscaling, options, inform[, match])

Find a matching-based symmetric scaling using the Hungarian algorithm.

The scaled matrix is such that the entry of maximum absolute value in each row and column is \(1.0\).

Parameters:

• m [integer ,in] :: number of rows in \(A\).

• n [integer ,in] :: number of columns in \(A\).

• ptr (n+1) [integer ,in] :: columns pointers for \(A\) (see CSC format).

• row (ptr(n+1)-1) [integer ,in] :: row indices for \(A\) (see CSC format).

• val (ptr(n+1)-1) [real ,in] :: non-zero values for \(A\) (see CSC format).

• rscaling (m) [real ,out] :: returns row scaling found by routine.

• cscaling (n) [real ,out] :: returns column scaling found by routine.

• options [hungarian_options ,in] :: controls behaviour of routine.

• inform [hungarian_inform ,out] :: returns information on execution of routine.

Options:

match (n) [integer ,out] :: returns matching found by routine. Row i is matched to column match(i), or is unmatched if match(i)=0.

Data-types

type  hungarian_options

Used to specify options to the routines hungarian_scale_sym() and hungarian_scale_unsym(). Components are automatically given default values in the definition of the type.

Please refer to the method section for details on how these parameters are used.

Type fields:

• % scale_if_singular [logical ,default=.false.] :: specifies behaviour for structurally singular matrices. If true, a partial scaling corresponding to a maximum cardinality matching will be returned. If false, an identity scaling is returned with an error code.

Note: If options%scale_if_singular=true, the resulting scaling will only be maximal for the matched rows/columns, and extreme care shuold be taken to ensure its use is meaningful!

type  hungarian_inform

Used to return information about the execution of the algorithm.

Type fields:

• % flag [integer ] :: gives the exit status of the algorithm (see table below)

• % matched [integer ] :: number of rows and columns that have been matched.

• % stat [integer ] :: holds the Fortran stat parameter in the event of an allocation failure (set to 0 otherwise).

Note: The number matched gives the structural rank of the matrix.

inform%flag	Return status
0	Success.
+1	Warning: Matrix is structurally rank-deficient. Only returned if options%scale_if_singular=true.
-1	Error: Allocation failed. Fortran stat value is returned in inform%stat.
-2	Error: Matrix is structurally rank-deficient. Only returned if options%scale_if_singular=false. Scaling vector(s) will be set to the identity, and a maximum cardinality matching will be returned in match(:) (if present).

Example

The following code shows an example usage of hungarian_scale_sym().

! examples/Fortran/scaling/hungarian_sym.f90 - Example code for SPRAL_SCALING
program hungarian_scale_sym_example
   use spral_scaling
   use spral_matrix_util, only : print_matrix, &
                                 SPRAL_MATRIX_REAL_SYM_INDEF
   implicit none

   ! Derived types
   type (hungarian_options)  :: options
   type (hungarian_inform)   :: inform

   ! Parameters
   integer, parameter :: wp = kind(0.0d0)

   ! Matrix data
   integer :: n, ptr(6), row(8)
   real(wp) :: val(8)

   ! Other variables
   integer :: match(5), i, j
   real(wp) :: scaling(5)

   ! Data for symmetric matrix:
   ! ( 2  1         )
   ! ( 1  4  1    8 )
   ! (    1  3  2   )
   ! (       2      )
   ! (    8       2 )
   n = 5
   ptr(1:n+1)        = (/ 1,        3,             6,      8,8,   9 /)
   row(1:ptr(n+1)-1) = (/ 1,   2,   2,   3,   5,   3,   4,   5   /)
   val(1:ptr(n+1)-1) = (/ 2.0, 1.0, 4.0, 1.0, 8.0, 3.0, 2.0, 2.0 /)
   write(*, "(a)") "Initial matrix:"
   call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)

   ! Perform symmetric scaling
   call hungarian_scale_sym(n, ptr, row, val, scaling, options, inform, &
      match=match)
   if(inform%flag<0) then
      write(*, "(a, i5)") "hungarian_scale_sym() returned with error ", &
         inform%flag
      stop
   endif

   ! Print scaling and matching
   write(*,"(a,10i10)")    'Matching:', match(1:n)
   write(*,"(a,10es10.2)") 'Scaling: ', scaling(1:n)

   ! Calculate scaled matrix and print it
   do i = 1, n
      do j = ptr(i), ptr(i+1)-1
         val(j) = scaling(i) * val(j) * scaling(row(j))
      end do
   end do
   write(*, "(a)") "Scaled matrix:"
   call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)

end program hungarian_scale_sym_example

The above code produces the following output:

Initial matrix:
Real unsymmetric matrix, dimension 5x5 with 10 entries.
1:   2.0000E+00   5.0000E+00
2:   1.0000E+00   4.0000E+00                             7.0000E+00
3:                1.0000E+00                2.0000E+00
4:                             3.0000E+00
5:                8.0000E+00                             2.0000E+00
Matching:         1         5         4         3         2
Row Scaling:   5.22E-01  5.22E-01  5.22E-01  5.22E-01  5.22E-01
Col Scaling:   9.59E-01  2.40E-01  6.39E-01  9.59E-01  2.74E-01
Scaled matrix:
Real unsymmetric matrix, dimension 5x5 with 10 entries.
1:   1.0000E+00   6.2500E-01
2:   5.0000E-01   5.0000E-01                             1.0000E+00
3:                1.2500E-01                1.0000E+00
4:                             1.0000E+00
5:                1.0000E+00                             2.8571E-01

Method

This algorithm is the same as used by the HSL package MC64. A scaling is derived from dual variables found during the solution of the below maximum product optimization problem using the Hungarian algorithm.

\[\begin{split}\max_{\sigma} & \prod_{i=1}^m\prod_{j=1}^n |a_{ij}|\sigma_{ij} & \\ \mathrm{s.t.} & \sum_{i=1}^m\sigma_{ij} = 1, & \forall j=1,n \\ & \sum_{j=1}^n\sigma_{ij} = 1, & \forall i=1,m \\ & \sigma_{ij} \in \{0,1\}.\end{split}\]

The array \(\sigma\) gives a matching of rows to columns.

By using the transformation

\[w_{ij} = \log c_j - \log |a_{ij}|,\]

where \(c_j = \max_i |a_{ij}|\), the maximum product problem in \(a_{ij}\) is replaced by a minimum sum problem in \(w_{ij}\) where all entries are positive. By standard optimization theory, there exist dual variables \(u\) and \(v\) corresponding to the constraints that satisfy the first order optimality conditions

\[\begin{split}w_{ij} - u_i - v_j = 0, & \mbox{if } \sigma_{ij }=1, \\ w_{ij} - u_i - v_j \ge 0, & \mbox{if } \sigma_{ij }=0.\end{split}\]

To obtain a scaling we define scaling matrices \(D_r\) and \(D_c\) as

\[ \begin{align}\begin{aligned}d^r_i = e^{u_i},\\d^c_i = e^{v_i}.\end{aligned}\end{align} \]

If a symmetric scaling is required, we average these as

\[d_i = \sqrt{d^r_id^c_i}.\]

By the first order optimality conditions, these scaling matrices guarantee that

\[\begin{split}d^r_i|a_{ij}|d^c_j = 1, && \mbox{if } \sigma_{ij}=1, \\ d^r_i|a_{ij}|d^c_j \le 1, && \mbox{if } \sigma_{ij}=0.\end{split}\]

To solve the minimum sum problem, the Hungarian algorithm maintains an optimal matching on a subset of the rows and columns. It proceeds to grow this set by finding augmenting paths from an unmatched row to an unmatched column. The algorithm is guaranteed to find the optimal solution in a fixed number of steps, but can be very slow as it may need to explore the full matrix a number of times equal to the dimension of the matrix. To minimize the solution time, a warmstarting heuristic is used to construct an initial optimal subset matching.

Further details are given in the following paper:

[3]

I.S. Duff and J. Koster. (1997). The design and use of algorithms for permuting large entries to the diagonal of sparse matrices. SIAM J. Matrix Anal. Applics. 20(4), pp 889–901. [Journal] [Preprint]