Version 6.0.2
with earlier contributions by Bob Miller & Mark Stankus
NCAlgebra is distributed under the terms of the BSD License:
Copyright (c) 2023, J. William Helton and Mauricio C. de Oliveira
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name NCAlgebra nor the names of its contributors may be
used to endorse or promote products derived from this software without
specific prior written permission.

This work was partially supported by the Division of Mathematical Sciences of the National Science Foundation.
The program was written by the authors with major earlier contributions from:
Considerable recent help came from
Other contributors include:
The beginnings of the program come from eran@slac.
NCPolyMonomial
in lex order.NCFromDigits
.NCGB
has been completely deprecated. Loading
NCGB
loads NCAlgebra
and NCGBX
instead.NCAlgebra is now distributed as a paclet!
Changed cannonical representation of noncommutative expressions
to allow for powers to be present in
NonCommutativeMultiply
.
 THE MOST NOTABLE LIKELY CONSTRUCTION THAT IS AFFECTED BY THIS CHANGE IS THE APPLICATION OF RULES BASED ON PATTERN MATCHING, WHICH NOW NEED TO EXPLICITLY TAKE INTO ACCOUNT THE PRESENCE OF EXPONENTS. SEE NOTES ON MANUAL FOR DETAILS ON HOW TO MITIGATE THE IMPACT OF THIS CHANGE. ALL NCALGEBRA COMMANDS HAVE BEEN REWRITTEN TO ACCOMODATE FOR THIS CHANGE IN REPRESENTATION.
NCTEST
renamed NCCORETEST
Tests must now be run as contexts:
e.g. << NCCORETEST`
instead of
<< NCCORETEST
NonCommutativeMultiply
: new functions NCExpandExponents and NCToList.
NCReplace
: new functions NCReplacePowerRule, NCExpandReplaceRepeated; NCExpandReplaceRepeatedSymmetric;
NCExpandReplaceRepeatedSelfAdjoint;
new option ApplyPowerRule
.
NCGBX
: NCMakeGB option
ReduceBasis
now defaults to True
.
NCCollect
: new function NCCollectExponents.
MatrixDecompositions
: functions GetLDUMatrices and GetFullLDUMatrices now produces low rank
matrices.
NCPoly
: new function NCPolyFromGramMatrixFactors.
NCPolyFullReduce
renamed NCPolyReduceRepeated.
NCPolyInterface
: new functions NCToRule, NCReduce, NCReduceRepeated, NCRationalToNCPoly, NCMonomialOrder, and NCMonomialOrderQ.
New utility functions SetCommutativeFunction, SetNonCommutativeFunction NCSymbolOrSubscriptExtendedQ, and NCNonCommutativeSymbolOrSubscriptExtendedQ.
The old C++
version of NCGB
is no
longer compatible with NCAlgebra
version 6.
Consider using NCGBX
instead.
No longer loads the package Notation
by default.
Controlled by the new option UseNotation
in NCOptions
.
Streamlined rules for NCSimplifyRational.
NC
and NCAlgebra
are now ContextsNCCollect
and NCStrongCollect
can handle
commutative variables.SetNonCommutativeHold
with
HoldAll
attribute can be used to set Symbols that have been
previously assigned values.NCWebInstall
and
NCWebUpdate
.NCGB
features are not fully supported yet, most notably NCProcess.CommuteEverything
(see important
notes in CommuteEverything).Transform
, Substitute
,
SubstituteSymmetric
, etc, have been replaced by the much
more reliable commands in the new package NCReplace.MatMult
has been replaced by NCDot. Alias MM
has been deprecated.x^3
expanding to x**x**x
, x^-1
expanding to
inv[x]
.x^T
expands to tp[x]
and x^*
expands to aj[x]
. Symbol T
is now
protected.This User Guide attempts to document the many improvements
introduced in NCAlgebra
Version 6. Please
be patient, as we move to incorporate the many recent changes into this
document.
See Reference Manual for a detailed description of the available commands.
There are also notebooks in the NC/DEMOS
directory that
accompany each of the chapters of this user guide.
Starting with Version 6, it is recommended that NCAlgebra be installed using our paclet distribution. Just type:
PacletInstall["https://github.com/NCAlgebra/NC/blob/master/NCAlgebra-6.0.2.paclet?raw=true"];
for the latest version.
In the near future we plan to submit paclets to the Wolfram paclet repository for easier updates.
Alternatively, you can download and install NCAlgebra as outlined in the section Manual Installation.
In Mathematica (notebook or text interface), type
<< NCAlgebra`
to load NCAlgebra
.
Advanced options for controlling the loading of NC
and
NCAlgebra
can be found in here and
here.
If you performed a manual installation of NCAlgebra you will need to type
<< NC`
before loading NCAlgebra.
If this step fails, your installation has problems (check out installation instructions in Manual Installation). If your installation is succesful, you will see a message like:
NC::Directory: You are using the version of NCAlgebra which is found in: "/your_home_directory/NC".
In the paclet version, it is no longer necessary to load the context
NC
before running NCAlgebra.Loading the context
NC
in the paclet version is however still supported for backward compatibility. It does nothing more than posting the message:NC::Directory: You are using a paclet version of NCAlgebra.
Extensive documentation is found at
and in the distribution directory
https://github.com/NCAlgebra/NC/DOCUMENTATION
which includes this document.
You may want to try some of the several demo files in the directory
DEMOS
after installing NCAlgebra
.
You can also run some tests to see if things are working fine.
There are 3 test sets which you can use to troubleshoot parts of NCAlgebra. The most comprehensive test set is run by typing:
<< NCCORETEST`
This will test the core functionality of NCAlgebra.
You can test functionality related to the package NCPoly
, including the new
NCGBX
package NCGBX
, by typing:
<< NCPOLYTEST`
Finally our Semidefinite Programming Solver NCSDP
can be tested with
<< NCSDPTEST`
We recommend that you restart the kernel before and after running tests. Each test takes a few minutes to run.
You can also call
<< NCPOLYTESTGB`
to perform extensive and long testing of NCGBX
.
If you performed a manual installation of NCAlgebra you will need to type
<< NC`
before before running any of the tests above.
Starting with Version 6, the old C++
version of our Groebner Basis Algorithm is no longer included. Consider
using NCGBX
.
This chapter provides a gentle introduction to some of the commands
available in NCAlgebra
.
If you want a living analog of this chapter just run the notebook
NC/DEMOS/1_MostBasicCommands.nb
.
Before you can use NCAlgebra
you first load it with the
following commands:
<< NC`
<< NCAlgebra`
If you installed the paclet version of NCAlgebra it is not necessary to load the context
NC
before loading otherNCAlgebra
packages. A dummy packageNC
is provided in case you would like to keep your NCAlgebra work compatible with previous versions.
In NCAlgebra
, the operator **
denotes
noncommutative multiplication. At present, single-letter lower
case variables are noncommutative by default and all others are
commutative by default. For example:
a**b-b**a
results in
a**b-b**a
while
A**B-B**A
A**b-b**A
both result in 0
.
One of Bill’s favorite commands is CommuteEverything
,
which temporarily makes all noncommutative symbols appearing in a given
expression to behave as if they were commutative and returns the
resulting commutative expression. For example:
CommuteEverything[a**b-b**a]
results in 0
. The command
EndCommuteEverything[]
restores the original noncommutative behavior.
One can make any symbol behave as noncommutative using
SetNonCommutative
. For example:
SetNonCommutative[A,B]
A**B-B**A
results in:
A**B-B**A
Likewise, symbols can be made commutative using
SetCommutative
. For example:
SetNonCommutative[A]
SetCommutative[B]
A**B-B**A
results in 0
. SNC
is an alias for
SetNonCommutative
. So, SNC
can be typed rather
than the longer SetNonCommutative
:
SNC[A];
A**a-a**A
results in:
-a**A+A**a
One can check whether a given symbol is commutative or not using
CommutativeQ
or NonCommutativeQ
. For
example:
CommutativeQ[B]
NonCommutativeQ[a]
both return True
.
WARNING: Prior to Version 6, noncommutative monomials would be stored in expanded form, without exponents. For example, the monomial
a**b**a^2**b
would be stored as
NonCommutativeMultiply[a, b, a, a, b]
The automatic expansion of powers of noncommutative symbols required overloading the behavior of the built in
Power
operator, which was interfering and causing much trouble when commutative polynomial operations were performed inside anNCAlgebra
session.Starting with Version 6, noncommutative monomials are represented with exponents. For instance, the same monomial above is now represented as
NonCommutativeMultiply[a, b, Power[a, 2], b]
Even if you type
a**b**a**a**b
, the repeated symbols get compressed to the compact representation with exponents. Exponents are now also used to represent the noncommutative inverse. See the notes in the next section.
The multiplicative identity is denoted Id
in the
program. At the present time, Id
is set to 1.
A symbol a
may have an inverse, which will be denoted by
inv[a]
. inv
operates as expected in most
cases.
For example:
inv[a]**a
inv[a**b]**a**b
both lead to Id = 1
and
a**b**inv[b]
results in a
.
WARNING: Starting with Version 6, the
inv
operator acts mostly as a wrapper forPower
. For example,inv[a]
internally evaluates toPower[a, -1]
. However,inv
is still available and used to display the noncommutative inverse of noncommutative expressions outside of a notebook environment. Beware thatinv[a**a]
now evaluates toPower[a, -2]
, hence certain patterns may no longer work. For example:NCReplace[inv[a**a], a**a -> b]
would produce
inv[b]
in previous versions ofNCAlgebra
but will fail in Version 6.
WARNING: Since Version 5,
inv
no longer automatically distributes over noncommutative products. If this more aggressive behavior is desired useSetOptions[inv, Distribute -> True]
. For exampleSetOptions[inv, Distribute -> True] inv[a**b]
returns
inv[b]**inv[a]
. ConverselySetOptions[inv, Distribute -> False] inv[a**b]
returns
inv[a**b]
.
tp[x]
denotes the transpose of symbol x
and
aj[x]
denotes the adjoint of symbol x
. Like
inv
, the properties of transposes and adjoints that
everyone uses constantly are built-in. For example:
tp[a**b]
leads to
tp[b]**tp[a]
and
tp[a+b]
returns
tp[a]+tp[b]
Likewise, tp[tp[a]] == a
and tp
for
anything for which CommutativeQ
returns True
is simply the identity. For example tp[5] == 5
,
tp[2 + 3I] == 2 + 3 I
, and tp[B] == B
.
Similar properties hold to aj
. Moreover
aj[tp[a]]
tp[aj[a]]
return co[a]
where co
stands for
complex-conjugate.
WARNING: Since Version 5 transposes (
tp
), adjoints (aj
), complex conjugates (co
), and inverses (inv
) in a notebook environment render as \(x^T\), \(x^*\), \(\bar{x}\), and \(x^{-1}\).tp
andaj
can also be input directly asx^T
andx^*
. For this reason the symbolT
is now protected inNCAlgebra
.
A trace like operator, tr
, was introduced in
v5.0.6. It is a commutative operator keeps its list of
arguments cyclicly sorted so that tr[b**a]
evaluates to
tr[a**b]
and that automatically distribute over sums so
that an expression like
tr[a**b - b**a]
always simplifies to zero. Also b**a**tr[b**a]
simplifies to
tr[a**b] a**b
because tr
is a commutative function. See SetCommutativeFunction.
A more interesting example is
expr = (a**b - b**a)^3
for which
NCExpand[tr[expr]]
evaluates to
3 tr[a^2**b^2**a**b] - 3 tr[a^2**b**a**b^2]
Use NCMatrixExpand to expand
tr
over matrices with noncommutative entries. For
example,
NCMatrixExpand[tr[{{a,b},{c,d}}]]
evaluates to
tr[a] + tr[d]
A key feature of symbolic computation is the ability to perform
substitutions. The Mathematica substitute commands,
e.g. ReplaceAll
(/.
) and
ReplaceRepeated
(//.
), are not reliable in
NCAlgebra
, so you must use our NC
versions of
these commands. For example:
NCReplaceAll[x**a**b,a**b->c]
results in
x**c
and
NCReplaceAll[tp[b**a]+b**a,b**a->c]
results in
c+tp[a]**tp[b]
Use NCMakeRuleSymmetric and NCMakeRuleSelfAdjoint to automatically create symmetric and self adjoint versions of your rules:
NCReplaceAll[tp[b**a]+b**a, NCMakeRuleSymmetric[b**a -> c]]
returns
c + tp[c]
WARNING: The change in internal representation introduced in Version 6, in which repeated letters in monomials are represented as powers, presents a new challenge to pattern matching for
NCAlgebra
expressions. For example, the seemingly innocent substitutionNCReplaceAll[a**b**b, a**b -> c, ApplyPowerRule -> False]
which in previous versions returned
c**b
will fail in Version 6. The reason for the failure is thata**b**b
is now internally represented asa**Power[b, 2]
, which does not matcha**b
. In order to make rules with exponents work in Version 6, they to be first modified by the new command NCReplacePowerRule as inNCReplaceAll[a**b**b, NCReplacePowerRule[a**b -> c], ApplyPowerRule -> False]
For convenience, when the option
ApplyPowerRule
is set toTrue
,NCReplacePowerRule
gets automatically applied by allNCReplace
family of functions. In this way,NCReplaceAll[a**b**b, a**b -> c] NCReplaceAll[a**b**b, a**b -> c, ApplyPowerRule -> True] NCReplaceAll[a**b**b, NCReplacePowerRule[a**b -> c], ApplyPowerRule -> False]
all return the more familiar result
c**b
in Version 6.
WARNING: NCReplacePowerRule and the option
ApplyPowerRule
may not work with the most exoteric replacements. See, for example, the note in section Inverses.
The difference between NCReplaceAll
and
NCReplaceRepeated
can be understood in the example:
NCReplaceAll[a**b^2, a**b -> a]
that results in
a**b
and
NCReplaceRepeated[a**b^2, a**b -> a]
that results in
a
Beside NCReplaceAll
and NCReplaceRepeated
we offer NCReplace
and NCReplaceList
, which
are analogous to the standard ReplaceAll
(/.
),
ReplaceRepeated
(//.
), Replace
and ReplaceList
. Note that one rarely uses
NCReplace
and NCReplaceList
.
With Version 6 we introduced NCExpandReplaceRepeated
,
NCExpandReplaceRepeatedSymmetric
,
and NCExpandReplaceRepeatedSelfAdjoint
,
which automate the tedious process of successive substitutions and
expansions that may be required to fully simplify expressions. For
example, consider the expression
expr = a**b^2-b^2**a
and the rule
rule = a**b -> a - b
for which
NCReplaceRepeated[expr, rule]
results in the expression
(a - b)**b - b^2**a
Note the presence of parenthesized terms resulting from the rule substitution. It is clear that after expanding
NCExpand[NCReplaceRepeated[expr, rule]]
to produce
a**b - b^2 - b^2**a
there are still terms that could be affected by the original replacement rule, that, if replaced again,
NCExpand[NCReplaceRepeated[expr, rule]]
NCExpand[NCReplaceRepeated[%, rule]]
would ultimately lead to a simpler expression
a - b - b^2 - b^2**a
This successive expansion and substitution process is automated in
NCExpandReplaceRepeated[expr, rule]
which produces the final expression
a - b - b^2 - b^2**a
in one shot.
See also the Sections Polynomials and
Rules and Advanced Rules and
Replacement for a deeper discussion on some issues involved with
rules and replacements in NCAlgebra
.
WARNING: The commands
Substitute
andTransform
have been deprecated in Version 5 in favor of the above nc versions ofReplace
.
The command NCExpand
expands noncommutative products.
For example:
NCExpand[(a+b)**x]
returns
a**x+b**x
Conversely, one can collect noncommutative terms involving same
powers of a symbol using NCCollect
. For example:
NCCollect[a**x+b**x,x]
recovers
(a+b)**x
NCCollect
groups terms by degree before collecting and
accepts more than one variable. For example:
expr = a**x+b**x+y**c+y**d+a**x**y+b**x**y
NCCollect[expr, {x}]
returns
y**c+y**d+(a+b)**x**(1+y)
and
NCCollect[expr, {x, y}]
returns
(a+b)**x+y**(c+d)+(a+b)**x**y
Note that the last term has degree 2 in x
and
y
and therefore does not get collected with the first order
terms.
The list of variables accepts tp
, aj
and
inv
, and
NCCollect[tp[x]**a**x+tp[x]**b**x+z,{x,tp[x]}]
returns
z+tp[x]**(a+b)**x
Alternatively one could use
NCCollectSymmetric[tp[x]**a**x+tp[x]**b**x+z,{x}]
to obtain the same result. A similar command, NCCollectSelfAdjoint, works with self-adjoint variables.
There is also a stronger version of collect called
NCStrongCollect
. NCStrongCollect
does not
group terms by degree. For instance:
NCStrongCollect[expr, {x, y}]
produces
y**(c+d)+(a+b)**x**(1+y)
Keep in mind that NCStrongCollect
often collects
more than one would normally expect.
WARNING: In Version 6 the commands NCExpandExponents and NCCollectExponents expands and collects exponents of expressions in noncommutative monomials. For example
NCExpandExponents[a**(b**a)^2**b]
produces
a**b**a**b**a**b
and
NCCollectExponents[a**b**a**b**a**b]
produces
(a**b)^3
NCAlgebra
provides some commands for noncommutative
polynomial manipulation that are similar to the native Mathematica
(commutative) polynomial commands. Some of the next commands required
the loading of the package
<< NCPolyInterface`
which provides an interface between NCAlgebra
and the
low-level package NCPoly.
For example:
SetCommutative[A, B]
expr = B + A y**x**y - 2 x
vars = NCVariables[expr]
returns
{x,y}
and
NCCoefficientList[expr, vars]
NCMonomialList[expr, vars]
NCCoefficientRules[expr, vars]
returns
{B, -2, A}
{1, x, y**x**y}
{1 -> B, x -> -2, y**x**y -> A}
Also for testing
NCMonomialQ[expr]
will return False
and
NCPolynomialQ[expr]
will return True
.
Another useful command is NCTermsOfDegree
, which will
returns an expression with terms of a certain degree. For instance:
NCTermsOfDegree[x**y**x - x^2**y + x**w + z**w, {x,y}, {2,1}]
returns x**y**x - x^2**y
,
NCTermsOfDegree[x**y**x - x^2**y + x**w + z**w, {x,y}, {0,0}]
returns z**w
, and
NCTermsOfDegree[x**y**x - x^2**y + x**w + z**w, {x,y}, {0,1}]
returns 0
.
A similar command is NCTermsOfTotalDegree
, which works
just like NCTermsOfDegree
but considers the total degree in
all variables. For example:
For example,
NCTermsOfTotalDegree[x**y**x - x^2**y + x**w + z**w, {x,y}, 3]
returns x**y**x - x^2**y
, and
NCTermsOfTotalDegree[x**y**x - x^2**y + x**w + z**w, {x,y}, 2]
returns 0
.
The above commands are based on special packages for efficiently storing and calculating with nc polynomials. Those packages are
NCPolynomial
: which
handles polynomials with noncommutative coefficients, andNCPoly
: which handles
polynomials with real coefficients.For example:
1 + y**x**y - A x
is a polynomial with commutative coefficients in \(x\) and \(y\), whereas
a**y**b**x**c**y - A x**d
is a polynomial with nc coefficients in \(x\) and \(y\), where the letters \(a\), \(b\), \(c\), and \(d\), are the nc coefficients. Of course
1 + y**x**y - A x
is a polynomial with nc coefficients if one considers only \(x\) as the variable of interest.
In order to take full advantage of NCPoly
and NCPolynomial
one would need
to convert an expression into those special formats. See the
sections on polynomials with
commutative coefficients and polynomials with
noncommutative coefficients in the mode advanced commands chapter and the
chapter on Gröebner basis for more information.
Details can be found in the package documentaions NCPolyInterface, NCPoly, and NCPolynomial.
One of the great challenges of noncommutative symbolic algebra is the
simplification of rational nc expressions. NCAlgebra
provides various algorithms that can be used for simplification and
general manipulation of nc rationals.
One such function is NCSimplifyRational
, which attempts
to simplify noncommutative rationals using a predefined set of rules.
For example:
expr = 1+inv[d]**c**inv[S-a]**b-inv[d]**c**inv[S-a+b**inv[d]**c]**b \
-inv[d]**c**inv[S-a+b**inv[d]**c]**b**inv[d]**c**inv[S-a]**b
NCSimplifyRational[expr]
leads to 1
. Of course the great challenge here is to
reveal well known identities that can lead to simplification. For
example, the two expressions:
expr1 = a**inv[1+b**a]
expr2 = inv[1+a**b]**a
and one can use NCSimplifyRational
to test such
equivalence by evaluating
NCSimplifyRational[expr1 - expr2]
which results in 0
or
NCSimplifyRational[expr1**inv[expr2]]
which results in 1
. NCSimplifyRational
works by transforming nc rationals. For example, one can verify that
NCSimplifyRational[expr2] == expr1
NCAlgebra
has a number of packages that can be used to
manipulate rational nc expressions. The packages:
NCGBX
perform calculations
with nc rationals using Gröbner basis, andNCRational
creates
state-space representations of nc rationals. This package is still
experimental.The package NCDiff
provide
functions for calculating derivatives and integrals of nc polynomials
and nc rationals.
<< NCDiff`
The main command is NCDirectionalD
which calculates
directional derivatives in one or many variables. For example, if:
expr = a**inv[1+x]**b + x**c**x
then
NCDirectionalD[expr, {x,h}]
returns
h**c**x + x**c**h - a**inv[1+x]**h**inv[1+x]**b
In the case of more than one variables
NCDirectionalD[expr, {x,h}, {y,k}]
takes the directional
derivative of expr
with respect to x
in the
direction h
and with respect to y
in the
direction k
. For example, if:
expr = x**q**x - y**x
then
NCDirectionalD[expr, {x,h}, {y,k}]
returns
h**q**x + x**q*h - y**h - k**x
A further example, if:
expr = x**a**x**b + x**c**x**d
then its directional derivative in the direction h
is
NCDirectionalD[expr, {x,h}]
which returns
h**a**x**b + x**a**h**b + h**c**x**d + x**c**h**d
The command NCGrad
calculates nc gradients1.
For example:
NCGrad[expr, x]
returns the nc gradient
a**x**b + b**x**a + c**x**d + d**x**c
A further example, if:
expr = x**a**x**b + x**c**y**d
is a function on variables x
and y
then
NCGrad[expr, x, y]
returns the nc gradient list
{a**x**b + b**x**a + c**y**d, d**x**c}
Version 5 introduced experimental support for
integration of nc polynomials. See NCIntegrate
.
NCAlgebra
has many commands for manipulating matrices
with noncommutative entries. Think block-matrices. Matrices are
represented in Mathematica using lists of lists. For
example
m = {{a, b}, {c, d}}
is a representation for the matrix
\[\left[\begin{array}{cc} a & b \\ c & d \end{array}\right]\]
The Mathematica command MatrixForm
output pretty
matrices. MatrixForm[m]
prints m
in a form
similar to the above matrix. Beware when copying and pasting parts of an
expression rendered byMatrixForm
because it may not execute
correctly. If in doubt, use FullForm
to reveal the contents
of the expression.
The experienced matrix analyst should always remember that the Mathematica convention for handling vectors is tricky.
{{1, 2, 4}}
is a 1x3 matrix or a row
vector;{{1}, {2}, {4}}
is a 3x1 matrix or a
column vector;{1, 2, 4}
is a vector but not
a matrix. Indeed whether it is a row or column vector depends
on the context. We advise not to use vectors.A useful command is NCInverse
,
which is akin to Mathematica’s Inverse
command and produces
a block-matrix inverse formula2 for an nc matrix. For
example
NCInverse[m]
returns
{{inv[a]**(1 + b**inv[d - c**inv[a]**b]**c**inv[a]), -inv[a]**b**inv[d - c**inv[a]**b]},
{-inv[d - c**inv[a]**b]**c**inv[a], inv[d - c**inv[a]**b]}}
or, using MatrixForm
,
NCInverse[m] // MatrixForm
returns
\[\begin{bmatrix} a^{-1} (1 + b (d - c a^{-1} b)^{-1} c a^{-1}) & -a^{-1} b (d - c a^{-1} b)^{-1} \\ -(d - c a^{-1} b)^{-1} c a^{-1} & (d - c a^{-1} b)^{-1} \end{bmatrix}\]
Note that a
and d - c**inv[a]**b
were
assumed invertible during the calculation.
Similarly, one can multiply matrices using NCDot
, which is similar to Mathematica’s
Dot
. For example
m1 = {{a, b}, {c, d}}
m2 = {{d, 2}, {e, 3}}
NCDot[m1, m2]
result in
{{a**d + b**e, 2 a + 3 b}, {c**d + d**e, 2 c + 3 d}}
Note that products of nc symbols appearing in the matrices are
multiplied using **
. Compare that with the standard
Dot
(.
) operator.
WARNING:
NCDot
replacesMatMult
, which is still available for backward compatibility but will be deprecated in future releases.
There are many new improvements with Version 5. For
instance, operators tp
, aj
, and
co
now operate directly over matrices. That is
aj[{{a,tp[b]},{co[c],aj[d]}}]
returns
{{aj[a],tp[c]},{co[b],d}}
In previous versions one had to use the special commands
tpMat
, ajMat
, and coMat
. Those
are still supported for backward compatibility.
See advanced matrix commands for
other useful matrix manipulation routines, such as NCMatrixExpand
, NCMatrixReplaceAll
, NCMatrixReplaceRepeated
,
etc, that allow one to work with matrices with symbolic noncommutative
entries.
Behind NCInverse
there are a host of linear algebra
algorithms which are available in the package:
NCMatrixDecompositions
:
implements versions of the \(LU\)
Decomposition with partial and complete pivoting, as well as \(LDL\) Decomposition which are suitable for
calculations with nc matrices. Those functions are based on the
templated algorithms from the package MatrixDecompositions
.For instance the function NCLUDecompositionWithPartialPivoting
can be used as
m = {{a, b}, {c, d}}
{lu, p} = NCLUDecompositionWithPartialPivoting[m]
which returns
lu = {{a, b}, {c**inv[a], d - c**inv[a]**b}}
p = {1, 2}
Using MatrixForm
:
MatrixForm[lu]
results in
\[\begin{bmatrix} a & b \\ c a^{-1} & d - c a^{-1} b \end{bmatrix}\]
The list p
encodes the sequence of permutations
calculated during the execution of the algorithm. The matrix
lu
contains the factors \(L\) and \(U\) in the way most common to numerical
analysts. These factors can be recovered using
{ll, uu} = GetFullLUMatrices[lu]
resulting in this case in
ll = {{1, 0}, {c**inv[a], 1}}
uu = {{a, b}, {0, d - c**inv[a]**b}}
Using MatrixForm
:
MatrixForm[ll]
MatrixForm[uu]
results in
\[L = \begin{bmatrix} 1 & 0 \\ c a^{-1} & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} a & b \\ 0 & d - c a^{-1} b \end{bmatrix}\]
To verify that \(M = L U\), input
m - NCDot[ll, uu]
which should return a zero matrix.
Note: if you are looking for efficiency, the
function GetLUMatrices
(also
GetLDUMatrices
) returns the
factors l
and u
as
SparseArrays
.
The default pivoting strategy prioritizes pivoting on simpler expressions. For instance,
m = {{a, b}, {1, d}}
{lu, p} = NCLUDecompositionWithPartialPivoting[m]
{ll, uu} = GetFullLUMatrices[lu]
result in the factors
ll = {{1, 0}, {a, 1}}
uu = {{1, d}, {0, b - a**d}}
and a permutation list
p = {2, 1}
which indicates that the number 1
, appearing in the
second row, was used as the pivot rather than the symbol a
appearing on the first row. Because of the permutation, to verify that
\(P M = L U\), input
m[[p]] - NCDot[ll, uu]
which should return a zero matrix. Note that in the above example the
permutation matrix \(P\) is never
constructed. Instead, the rows of \(M\)
are directly permuted using Mathematica’s Part
([[]]
) command. Of course, if one prefers to work with
permutation matrices, they can be easily obtained by permuting the rows
of the identity matrix as in the following example
p = {2, 1, 3}
IdentityMatrix[3][[p]] // MatrixForm
to produce
\[\begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}\]
Likewise
m = {{a + b, b}, {c, d}}
{lu, p} = NCLUDecompositionWithPartialPivoting[m]
{ll, uu} = GetFullLUMatrices[lu]
returns
p = {2, 1}
ll = {{1, 0}, {(a + b)**inv[c], 1}}
uu = {{c, d}, {0, b - (a + b)**inv[c]**d}}
showing that the simpler expression c
was taken
as a pivot instead of a + b
.
The function NCLUDecompositionWithPartialPivoting
is the
one that is used by NCInverse
.
Another factorization algorithm is NCLUDecompositionWithCompletePivoting
,
which can be used to calculate the symbolic rank of nc matrices. For
example
m = {{2 a, 2 b}, {a, b}}
{lu, p, q, rank} = NCLUDecompositionWithCompletePivoting[m]
returns the left and right permutation lists
p = {2, 1}
q = {1, 2}
and rank
equal to 1
. Note that
p = {2, 1}
and q = {1,2}
tell us that the
element that was pivoted on was the symbol a
, which is the
first entry of the second row, rather then 2 a
, which is
the first entry of the first row, because a
is
simpler than 2 a
. The \(L\) and \(U\) factors can be obtained as before
using
{ll, uu} = GetFullLUMatrices[lu]
to get
ll = {{1, 0}, {2, 1}}
uu = {{a, b}, {0, 0}}
Using MatrixForm
:
MatrixForm[ll]
MatrixForm[uu]
results in
\[L = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} a & b \\ 0 & 0 \end{bmatrix}\]
In this case, to verify that \(P M Q = L U\) input
NCDot[ll, uu] - m[[p, q]]
which should return a zero matrix. As with partial pivoting, the
permutation matrices \(P\) and \(Q\) are never constructed. Instead we used
Part
([[]]
) to permute both rows and
columns.
Finally NCLDLDecomposition
computes
the \(LDL^T\) decomposition of
symmetric symbolic nc matrices. For example
m = {{a, b}, {b, c}}
{ldl, p, s, rank} = NCLDLDecomposition[m]
returns ldl
, which contain the factors, and
p = {1, 2}
s = {1, 1}
rank = 2
The list p
encodes left and right permutations,
s
is a list specifying the size of the diagonal blocks
(entries can be either 1 or 2). The factors can be obtained using GetLDUMatrices
as in
{ll, dd, uu} = GetFullLDUMatrices[ldl, s]
which in this case returns
ll = {{1, 0}, {b**inv[a], 1}}
dd = {{a, 0}, {0, c - b**inv[a]**b}}
uu = {{1, inv[a]**b}, {0, 1}}}
Because \(P M P^T = L D L^T\),
NCDot[ll, dd, uu] - m[[p, p]]
is the zero matrix and \(U = L^T\).
NCLDLDecomposition
works only on symmetric matrices and,
whenever possible, will make invertibility and symmetry assumptions on
variables so that it can run successfully. If not possible it will warn
the users.
WARNING: Prior versions contained the command
NCLDUDecomposition
which was deprecated in Version 5 as its functionality is now provided byNCLDLDecomposition
, with a slightly different syntax.
NCMatrixReplaceAll
and
NCMatrixReplaceRepeated
are special versions of NCReplaceAll
and NCReplaceRepeated
that take
extra steps to preserve matrix consistency when replacing expressions
with nc matrices. For example, with
m1 = {{a, b}, {c, d}}
m2 = {{d, 2}, {e, 3}}
and
M = {{a11,a12}}
the call
MM = NCMatrixReplaceRepeated[M, {a11 -> m1, a12 -> m2}]
produces as a result the matrix
{{a, b, d, 2}, {c, d, e, 3}}
or, using MatrixForm
:
MatrixForm[MM]
to obtain
\[\begin{bmatrix} a & b & d & 2 \\ c & d & e & 3 \end{bmatrix}\]
Note how the symbols were treated as block-matrices during the substitution. As a second example, with
M = {{a11, 0}, {0, a22}}
the command
MM = NCMatrixReplaceRepeated[M, {a11 -> m1, a22 -> m2}]
produces the matrix
{{a, b, 0, 0}, {c, d, 0, 0}, {0, 0, d, 2}, {0, 0, e, 3}}
or, using MatrixForm
:
MatrixForm[MM]
to obtain
\[\begin{bmatrix} a & b & 0 & 0 \\ c & d & 0 & 0 \\ 0 & 0 & d & 2 \\ 0 & 0 & e & 3 \end{bmatrix}\]
in which the 0
blocks were automatically expanded to fit
the adjacent block matrices.
Another feature of NCMatrixReplace
and its variants is
its ability to withhold evaluation until all matrix substitutions have
taken place. For example,
NCMatrixReplaceAll[x**y + y, {x -> m1, y -> m2}]
produces
{{d + a**d + b**e, 2 + 2 a + 3 b},
{e + c**d + d**e, 3 + 2 c + 3 d}}
Finally, NCMatrixReplace
substitutes
NCInverse
for inv
so that, for instance, the
result of
rule = {x -> m1, y -> m2, id -> IdentityMatrix[2], z -> {{id,x},{x,id}}}
NCMatrixReplaceRepeated[inv[z], rule]
coincides with
NCInverse[ArrayFlatten[{{IdentityMatrix[2], m1}, {m1, IdentityMatrix[2]}}]]
The closest related demo to the material in this section is
NC/DEMOS/NCConvexity.nb
.
When working with nc quadratics it is useful to be able to “factor” the quadratic into the following form
\[ q(x) = c + s(x) + l(x) M r(x) \]
where \(s\) is linear \(x\) and \(l\) and \(r\) are vectors and \(M\) is a matrix. Load the package
<< NCQuadratic`
and use the command NCToNCQuadratic
to factor an nc
polynomial into the the above form:
vars = {x, y};
expr = tp[x]**a**x**d + tp[x]**b**y + tp[y]**c**y + tp[y]**tp[b]**x**d;
{const, lin, left, middle, right} = NCToNCQuadratic[expr, vars];
which returns
left = {tp[x],tp[y]}
right = {y, x**d}
middle = {{a,b}, {tp[b],c}}
and zero const
and lin
. The format for the
linear part lin
will be discussed later in Section Linear. Note that coefficients of an nc quadratic may
also appear on the left and right vectors, as d
did in the
above example. Conversely, NCQuadraticToNC
converts a list
with factors back to an nc expression as in:
NCQuadraticToNC[{const, lin, left, middle, right}]
which results in
(tp[x]**b + tp[y]**c)**y + (tp[x]**a + tp[y]**tp[b])**x**d
An interesting application is the verification of the domain in which an nc rational function is convex. This uses the second directional derivative, called the Hessian. Take for example the quartic
expr = x^4;
and calculate its noncommutative directional Hessian
hes = NCHessian[expr, {x, h}]
This command returns
2 h^2**x^2 + 2 h**x**h**x + 2 h**x^2**h + 2 x**h^2**x + 2 x**h**x**h + 2 x^2**h^2
which is quadratic in the direction h
. The decomposition
of the nc Hessian using NCToNCQuadratic
{const, lin, left, middle, right} = NCToNCQuadratic[hes, {h}];
produces
left = {h, x**h, x^2**h}
right = {h**x^2, h**x, h}
middle = {{2, 2 x, 2 x^2},{0, 2, 2 x},{0, 0, 2}}
Note that the middle matrix
\[\begin{bmatrix} 2 & 2 x & 2 x^2 \\ 0 & 2 & 2 x \\ 0 & 0 & 2 \end{bmatrix}\]
is not symmetric, as one might have expected. The command NCQuadraticMakeSymmetric
can fix that and produce a symmetric decomposition. For the above
example
{const, lin, sleft, smiddle, sright} =
NCQuadraticMakeSymmetric[{const, lin, left, middle, right},
SymmetricVariables -> {x, h}]
results in
sleft = {x^2**h, x**h, h}
sright = {h**x^2, h**x, h}
middle = {{0, 0, 2}, {0, 2, 2 x}, {2, 2 x, 2 x^2}}
in which middle
is the symmetric matrix
\[\begin{bmatrix} 0 & 0 & 2 \\ 0 & 2 & 2 x \\ 2 & 2 x & 2 x^2 \end{bmatrix}\]
Note the argument SymmetricVariables -> {x,h}
which
tells NCQuadraticMakeSymmetric
to consider x
and y
as symmetric variables. Because the
middle
matrix is never positive semidefinite for any
possible value of \(x\) the
conclusion3 is that the nc quartic \(x^4\) is not convex.
The production of such symmetric quadratic decompositions is
automated by the convenience command NCMatrixOfQuadratic
. Verify
that
{sleft, smiddle, sright} = NCMatrixOfQuadratic[hes, {h}]
automatically assumes that both x
and h
are
symmetric variables and produces suitable left and right vectors as well
as a symmetric middle matrix. Now we illustrate the application of such
command to checking the convexity region of a noncommutative rational
function.
If one is interested in checking convexity of nc rationals the
package NCConvexity
has
functions that automate the whole process, including the calculation of
the Hessian and the middle matrix, followed by the diagonalization of
the middle matrix as produced by NCLDLDecomposition
.
For example, the commands evaluate the nc Hessian and calculates its quadratic decomposition
expr = (x + b**y)**inv[1 - a**x**a + b**y + y**b]**(x + y**b);
{left, middle, right} = NCMatrixOfQuadratic[NCHessian[expr, {x, h}], {h}];
The resulting middle matrix can be factored using
{ldl, p, s, rank} = NCLDLDecomposition[middle];
{ll, dd, uu} = GetLDUMatrices[ldl, s];
which produces the diagonal factors
\[\begin{bmatrix} 2 (1 + b y + y b - a x a)^{-1} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\]
which indicates the the original nc rational is convex whenever
\[ (1 + b y + y b - a x a)^{-1} \succeq 0 \]
or, equivalently, whenever
\[ 1 + b y + y b - a x a \succeq 0 \]
The above sequence of calculations is automated by the command NCConvexityRegion
as in
<< NCConvexity`
NCConvexityRegion[expr, {x}]
which results in
{2 inv[1 + b**y + y**b - a**x**a], 0}
which correspond to the diagonal entries of the LDL decomposition of the middle matrix of the nc Hessian.
In this chapter we describe some more advance features and commands.
If you want a living version of this chapter just run the notebook
NC/DEMOS/2_MoreAdvancedCommands.nb
.
<< NC`
<< NCAlgebra`
Substitution is a key feature of Symbolic computation. We will now
discuss some issues related to Mathematica’s implementation of rules and
replacements that can affect the behavior of NCAlgebra
expressions.
ReplaceAll
(/.
) and ReplaceRepeated
(//.
)
often failThe first issue is related to how Mathematica performs substitutions, which is through pattern matching. For a rule to be effective if has to match the structural representation of an expression. That representation might be different than one would normally think based on the usual properties of mathematical operators. For example, one would expect the rule:
rule = 1 + x_ -> x
to match all the expressions bellow:
1 + a
1 + 2 a
1 + a + b
1 + 2 a * b
so that
expr /. rule
with expr
taking the above expressions would result
in:
a
2 a
a + b
2 a * b
Indeed, Mathematica’s attribute Flat
does precisely
that. Note that this is still structural matching, not
mathematical matching, since the pattern 1 + x_
would not match an integer 2
, even though one could write
2 = 1 + 1
!
Unfortunately, **
, which is the
NonCommutativeMultiply
operator, is not
Flat
4. This is the reason why substitution
based on a simple rule such as:
rule = a**b -> c
so that
expr /. rule
will work for some expr
like
1 + 2 a**b /. rule
resulting in
1 + 2 c
but will fail to produce the expected result in cases like:
a**b**c /. rule
or
c**a**b /. rule
c**a**b**d /. rule
1 + 2 a**b**c /. rule
That’s what the NCAlgebra
family of replacement
functions discussed in the next section are made for.
NCReplace
Continuing with the example in the previous section, the calls
NCReplaceAll[a**b**c, rule]
NCReplaceAll[ c**a**b, rule ]
NCReplaceAll[c**a**b**d, rule]
NCReplaceAll[1 + 2 a**b**c, rule ]
produce the results one would expect:
c^2
c^2
c^2**d
1 + 2 c^2
For this reason, when substituting in NCAlgebra
it is
always safer to use functions from the NCReplace
package rather than
the corresponding Mathematica Replace
family of functions.
Unfortunately, this comes at a the expense of sacrificing the standard
operators /.
(ReplaceAll
) and //.
(ReplaceRepeated
), which cannot be safely overloaded,
forcing one to use the full names NCReplaceAll
and
NCReplaceRepeated
.
On the same vein, the following substitution rule
NCReplaceAll[2 a**b + c, 2 a -> b]
will return 2 a**b + c
intact since
FullForm[2 a**b]
is actually
Times[2, NonCommutativeMuliply[a, b]]
which is not structurally related to
FullForm[2 a]
which is
Times[2, a]
Of course, in this case a simple solution is to use the alternative rule:
NCReplaceAll[2 a**b + c, a -> b / 2]
which results in b^2 + c
, as one might expect.
Starting with Version 6, NCAlgebra
stores repeated symbols in noncommutative monomials using powers. This
means that
expr = a**a**a**b**a**a**b**a**b
FullForm[expr]
is internally stored as
NonCommutativeMultiply[a^3, b, a^2, b, a, b]
so that a replacement such as
NCReplaceAll[expr, a**b -> c, ApplyPowerRule -> False]
will result in
a^3**b**a^2**b**c
Note how the rule fails to match the a**b
in the terms
a**b^2
and a**b^3
. This situation might be
familiar to an experienced Mathematica user who is a aware of the
difference between structural pattern matching and mathematical
matching5.
As a convenience for NCAlgebra
users, Version
6 provides the function NCReplacePowerRule that modifies a user’s
rule in order to accomplish matching of symbols in
NCAlgebra
expressions including powers. For example, for
the same expr
above, the following replacement
NCReplaceAll[expr, NCReplacePowerRule[a**b -> c], ApplyPowerRule -> False]
will produce the more familiar
a^2**c**a^2**b**a**b
after matching a**b
in the term a^3**b
,
and
NCReplaceRepeated[expr, NCReplacePowerRule[a**b -> c], ApplyPowerRule -> False]
will produce
a^2**c**a**c^2
after matching a**b
in a^3**b
,
a^2**b
, and a**b
.
The command NCReplacePowerRule
works by modifying
noncommutative patterns with symbols appearing at the edges of monomials
to account for the potential presence of Power
in a
noncommutative monomial. For example,
NCReplacePowerRule[a**c**b -> d]
produces the modified rule6
a^n_.**c**b^m_. -> a^(n-1)**d**b^(m-1)
which can successfully match powers of the symbols a
and
b
appearing in the monomial a**c**b
.
The application of NCReplacePowerRule
can also be done
by invoking the option ApplyPowerRule
, which is available
in all functions of the package NCReplace. Currently,
ApplyPowerRule
is set to True
by default so
that the commands
NCReplaceRepeated[expr, a**b -> c]
NCReplaceRepeated[expr, a**b -> c, ApplyPowerRule -> True]
do the same thing as
NCReplaceRepeated[expr, NCReplacePowerRule[a**b -> c], ApplyPowerRule -> False]
This option can also be turned off globally for all calls to the
NCReplace
family of functions in a NCAlgebra
session by calling
SetOptions[NCReplace, ApplyPowerRule -> False];
After that, all calls to NCReplace
,
NCReplaceAll
, NCReplaceRepeated
, and related
function, will be done with the option
ApplyPowerRule -> False
automatically.
To revert to the default behavior just set
SetOptions[NCReplace, ApplyPowerRule -> True];
Block
and Module
A second more esoteric issue related to substitution in
NCAlgebra
does not have a clean solution. It is also one
that usually lurks in hidden pieces of code and can be very difficult to
spot. We have been victim of such “bugs” many times. Luckily it only
affect advanced users that are using NCAlgebra
inside their
own functions using Mathematica’s Block
and
Module
constructions. It is also not a real bug, since it
follows from some often not well understood issues with the usage of Block
versus Module
. Our goal here is therefore not to
fix the issue but simply to alert advanced users of its
existence. Start by first revisiting the following example from the Mathematica
documentation. Let
m = i^2
and run
Block[{i = a}, i + m]
which returns the “expected”
a + a^2
versus
Module[{i = a}, i + m]
which returns the “surprising”
a + i**i
The reason for this behavior is that Block
effectively
evaluates i
as a local variable and evaluates
m
using whatever values are available at the time of
evaluation, whereas Module
only evaluates the symbol
i
which appears explicitly in i + m
and not m
using the local value of i = a
. This
can lead to many surprises when using rules and substitution inside a
Module
. For example:
Block[{i = a}, i_ -> i]
will return
i_ -> a
whereas
Module[{i = a}, i_ -> i]
will return
i_ -> i
More devastating for NCAlgebra
is the fact that
Module
will hide local definitions from rules, which will
often lead to disaster if local variables need to be declared
noncommutative. Consider for example the trivial definitions for
F
and G
below:
F[exp_] := Module[
{rule, aa, bb},
SetNonCommutative[aa, bb];
rule = aa_**bb_ -> bb**aa;
NCReplaceAll[exp, rule]
]
G[exp_] := Block[
{rule, aa, bb},
SetNonCommutative[aa, bb];
rule = aa_**bb_ -> bb**aa;
NCReplaceAll[exp, rule]
]
Their only difference is that one is defined using a
Block
and the other is defined using a Module
.
The task is to apply a rule that flips the noncommutative
product of their arguments, say, x**y
, into
y**x
. The problem is that only one of those definitions
work “as expected.” Indeed, verify that
G[x**y]
returns the “expected”
y**x
whereas
F[x**y]
returns
x y
which completely destroys the noncommutative product. The reason for
the catastrophic failure of the definition of F
, which is
inside a Module
, is that the letters aa
and
bb
appearing in rule
are not treated as
the local symbols aa
and bb
7.
For this reason, the right-hand side of the rule rule
involves the global symbols aa
and bb
, which
are, in the absence of a declaration to the contrary, commutative. On
the other hand, the definition of G
inside a
Block
makes sure that aa
and bb
are evaluated with whatever value they might have locally at the time of
execution.
The above subtlety often manifests itself partially, sometimes
causing what might be perceived as some kind of erratic
behavior. Indeed, if one had used symbols that were already
declared globally noncommutative by NCAlgebra
, such as
single small cap roman letters as in the definition:
H[exp_] := Module[
{rule, a, b},
SetNonCommutative[a, b];
rule = a_**b_ -> b**a;
NCReplaceAll[exp, rule]
]
then calling
H[x**y]
works “as expected,” even if for the wrong reasons!
Another possible “fix” is to use a delayed rule, as in:
H[exp_] := Module[
{rule, aa, bb},
SetNonCommutative[aa, bb];
rule = aa_**bb_ :> bb**aa;
NCReplaceAll[exp, rule]
]
with which
H[x**y]
would also work because the evaluation of the right-hand side of the rule is delayed until the time of its application.
Starting at Version 5 the operators **
and inv
apply also to matrices. However, in order for
**
and inv
to continue to work as full fledged
operators, the result of multiplications or inverses of matrices is held
unevaluated until the user calls NCMatrixExpand
. This is in the
the same spirit as good old fashion commutative operations in
Mathematica.
For example, with
m1 = {{a, b}, {c, d}}
m2 = {{d, 2}, {e, 3}}
the call
m1**m2
results in
{{a, b}, {c, d}}**{{d, 2}, {e, 3}}
Upon calling
m1**m2 // NCMatrixExpand
the matrix product evaluation takes place returning
{{a**d + b**e, 2a + 3b}, {c**d + d**e, 2c + 3d}}
which is what would have arisen from calling
NCDot[m1,m2]
8. Likewise
inv[m1]
results in
inv[{{a, b}, {c, d}}]
and
inv[m1] // NCMatrixExpand
returns the evaluated result
{{inv[a]**(1 + b**inv[d - c**inv[a]**b]**c**inv[a]), -inv[a]**b**inv[d - c**inv[a]**b]},
{-inv[d - c**inv[a]**b]**c**inv[a], inv[d - c**inv[a]**b]}}
or, using MatrixForm
:
inv[m1] // NCMatrixExpand // MatrixForm
returns
\[\begin{bmatrix} a^{-1} (1 + b (d - c a^{-1} b)^{-1} c a^{-1}) & -a^{-1} b (d - c a^{-1} b)^{-1} \\ -(d - c a^{-1} b)^{-1} c a^{-1} & (d - c a^{-1} b)^{-1} \end{bmatrix}\]
A less trivial example is
m3 = m1**inv[IdentityMatrix[2] + m1] - inv[IdentityMatrix[2] + m1]**m1
that returns
-inv[{{1 + a, b}, {c, 1 + d}}]**{{a, b}, {c, d}} +
{{a, b}, {c, d}}**inv[{{1 + a, b}, {c, 1 + d}}]
Expanding
NCMatrixExpand[m3]
results in
{{b**inv[b - (1 + a)**inv[c]**(1 + d)] - inv[c]**(1 + (1 + d)**inv[b -
(1 + a)**inv[c]**(1 + d)]**(1 + a)**inv[c])**c - a**inv[c]**(1 + d)**inv[b -
(1 + a)**inv[c]**(1 + d)] + inv[c]**(1 + d)**inv[b - (1 + a)**inv[c]**(1 + d)]**a,
a**inv[c]**(1 + (1 + d)**inv[b - (1 + a)**inv[c]**(1 + d)]**(1 + a)**inv[c]) -
inv[c]**(1 + (1 + d)**inv[b - (1 + a)**inv[c]**(1 + d)]**(1 + a)**inv[c])**d -
b**inv[b - (1 + a)**inv[c]**(1 + d)]**(1 + a)**inv[c] + inv[c]**(1 + d)**inv[b -
(1 + a)**inv[c]**(1 + d)]** b},
{d**inv[b - (1 + a)**inv[c]**(1 + d)] - (1 + d)**inv[b - (1 + a)**inv[c]**(1 + d)] -
inv[b - (1 + a)**inv[c]**(1 + d)]**a + inv[b - (1 + a)**inv[c]**(1 + d)]**(1 + a),
1 - inv[b - (1 + a)**inv[c]**(1 + d)]**b - d**inv[b - (1 + a)**inv[c]**(1 + d)]**(1 +
a)**inv[c] + (1 + d)**inv[b - (1 + a)**inv[c]**(1 + d)]**(1 + a)**inv[c] +
inv[b - (1 + a)**inv[c]**(1 + d)]**(1 + a)**inv[c]**d}}
and, finally,
NCMatrixExpand[m3] // NCSimplifyRational
returns
{{0, 0}, {0, 0}}
as expected.
Plus
(+
) and matricesMathematica’s choice of treating lists and matrix indistinctively can
cause much trouble when mixing **
with Plus
(+
) operator. The reason for that is that Plus
is Listable
, and an expression like
z + m1
where m1
is an array and z
is not, is
automatically expanded into
{{a + z, b + z}, {c + z, d + z}}
or, using MatrixForm
:
z + m1 // MatrixForm
resulting in
\[\begin{bmatrix} a + z & b + z \\ c + z & d + z \end{bmatrix}\]
Because of this “feature”, the expression
m1**m2 + m2 // NCMatrixExpand
evaluates to the “wrong” result
{{{{d + a**d + b**e, 2 a + 3 b + d}, {d + c**d + d**e, 2 c + 4 d}},
{{2 + a**d + b**e, 2 + 2 a + 3 b}, {2 + c**d + d**e, 2 + 2 c + 3 d}}},
{{{e + a**d + b**e, 2 a + 3 b + e}, {e + c**d + d**e, 2 c + 3 d + e}},
{{3 + a**d + b**e, 3 + 2 a + 3 b}, {3 + c**d + d**e, 3 + 2 c + 3 d}}}}
which is different than the “correct” result
{{d + a**d + b**e, 2 + 2 a + 3 b},
{e + c**d + d**e, 3 + 2 c + 3 d}}
which is returned by either
NCMatrixExpand[m1**m2] + m2
or
NCDot[m1, m2] + m2
The reason for the behavior is that m1**m2
is
essentially treated as a scalar (it does not have head
List
) and therefore gets added entrywise to m2
before NCMatrixExpand
has a chance to evaluate the
**
product.
There are no easy fixes for this problem, which affects not only
NCAlgebra
but any similar type of matrix product evaluation
in Mathematica. With NCAlgebra
, a better option is to use
NCMatrixReplaceAll
or NCMatrixReplaceRepeated
.
As seen in the section Replace with
matrices, the command
NCMatrixReplaceAll[x**y + y, {x -> m1, y -> m2}]
produces the expected result:
{{d + a**d + b**e, 2 + 2 a + 3 b},
{e + c**d + d**e, 3 + 2 c + 3 d}}
Note that the above behavior might be erratic, since in an expression like
m1**m2 + m2**m1
in which **
is kept unevaluated as in
{{a, b}, {c, d}}**{{d, 2}, {e, 3}} + {{d, 2}, {e, 3}}**{{a, b}, {c, d}}
the command
m1**m2 + m2**m1 // NCMatrixExpand
expands to the expected result
{{2 c + a**d + b**e + d**a, 2 a + 3 b + 2 d + d**b},
{3 c + c**d + d**e + e**a, 2 c + 6 d + e**b}}
The package NCPoly
provides
an efficient structure for storing and operating with noncommutative
polynomials with commutative coefficients. There are two main goals:
Those two properties allow for an efficient implementation of
NCAlgebra
’s noncommutative Gröbner basis algorithm, new in
Version 5, without the use of auxiliary accelerating
C
code, as in NCGB
. See Noncommutative Gröbner Basis.
Before getting into details, to see how much more efficient
NCPoly
is when compared with standard
NCAlgebra
objects try
Table[Timing[NCExpand[(1 + x)^i]][[1]], {i, 0, 15, 5}]
which would typically return something like
{0.000088, 0.001074, 0.017322, 0.240704, 3.61492}
whereas the equivalent
<< NCPoly`
Table[Timing[(1 + NCPolyMonomial[{x}, {x}])^i][[1]], {i, 0, 15, 5}]
would return
{0.00097, 0.001653, 0.002208, 0.003908, 0.004306}
The best way to work with NCPoly
in
NCAlgebra
is by loading the package NCPolyInterface
:
<< NCPolyInterface`
which provides the commands NCToNCPoly
and NCPolyToNC
to convert nc expressions
back and forth between NCAlgebra
and
NCPoly
.
For example
vars = {x, y, z};
p = NCToNCPoly[1 + x**x - 2 x**y**z, vars]
converts the polynomial 1 + x**x - 2 x**y**z
from the
standard NCAlgebra
format into an NCPoly
object. The result in this case is the NCPoly
object
NCPoly[{1, 1, 1}, <|{0, 0, 0, 0} -> 1, {0, 0, 2, 0} -> 1, {1, 1, 1, 5} -> -2|>]
Conversely the command NCPolyToNC
converts an
NCPoly
back into NCAlgebra
format. For
example
NCPolyToNC[p, vars]
returns
1 + x^2 - 2 x**y**z
as expected. Note that an NCPoly
object does not store
symbols, but rather a representation of the polynomial based on
specially encoded monomials. This is the reason why one should provide
vars
as an argument to NCPolyToNC
.
Alternatively, one could construct the same NCPoly
object by calling NCPoly
directly as in
NCPoly[{1, 1, -2}, {{}, {x, x}, {x, y, z}}, vars]
In this syntax the first argument is a list of coefficients,
the second argument is a list of monomials, and the third is
the list of variables. Monomials are given as lists,
with {}
being equivalent to a constant 1
.
The particular coefficients in the NCPoly
object depend
not only on the polynomial being represented but also on the
ordering implied by the sequence of symbols in the list of
variables vars
. For example:
vars = {{x}, {y, z}};
p = NCToNCPoly[1 + x**x - 2 x**y**z, vars]
produces:
NCPoly[{1, 2}, <|{0, 0, 0} -> 1, {0, 2, 0} -> 1, {2, 1, 5} -> -2|>
The sequence of braces in the list of variables encodes the
ordering to be used for sorting NCPoly
s. Orderings
specify how monomials should be ordered, and is discussed in detail in
Noncommutative Gröbner Basis.
We provide the convenience commands NCMonomialOrder
, NCMonomialOrderQ
and NCPolyDisplayOrder
to help
building and visualizing orderings.
NCPolyDisplayOrder
prints the polynomial ordering implied by a list of symbols. For
example
NCPolyDisplayOrder[{x,y,z}]
prints out
\[x \ll y \ll z\]
and
NCPolyDisplayOrder[{{x},{y,z}}]
prints out
\[x \ll y < z\]
from where you can see that grouping variables inside braces induces a graded type ordering, as discussed in Noncommutative Gröbner Basis.
With NCMonomialOrder
you
can basically dispense with the “outer” braces. For example
NCMonomialOrder[{x},{y,z}]
returns
{{x},{y,z}}
It also validates the resulting ordering using NCMonomialOrderQ
so that
NCMonomialOrder[{x},{y,{z}}]
will fail and print an error message. Note that
NCMonomialOrder[x,y,z]
returns
{{x},{y},{z}}
which corresponds to
\[x \ll y \ll z\]
and
NCMonomialOrder[{x, y, z}]
returns
{{x,y,z}}
which corresponds to
\[x < y < z\]
There is also a special constructor for monomials. For example
vars = NCMonomialOrder[{x}, {y,z}];
NCPolyMonomial[{y,x}, vars]
NCPolyMonomial[{x,y}, vars]
return the monomials corresponding to \(y x\) and \(x y\) using the ordering \(x \ll y < z\).
Operations on NCPoly
objects result in another
NCPoly
object that is always expanded. For example:
vars = {{x}, {y, z}};
1 + NCPolyMonomial[{x, y}, vars] - 2 NCPolyMonomial[{y, x}, vars]
returns
NCPoly[{1, 2}, <|{0, 0, 0} -> 1, {1, 1, 1} -> 1, {1, 1, 3} -> -2|>]
and
p = (1 + NCPolyMonomial[{x}, vars]**NCPolyMonomial[{y}, vars])^2
returns
NCPoly[{1, 2}, <|{0, 0, 0} -> 1, {1, 1, 1} -> 2, {2, 2, 10} -> 1|>]
WARNING:
NCPoly
s constructed from different orderings cannot be combined or otherwise operated.
Another convenience function is NCPolyDisplay
which
returns a list with the monomials appearing in an NCPoly
object. For example:
NCPolyDisplay[p, vars]
returns
{x.y.x.y, 2 x.y, 1}
The reason for displaying an NCPoly
object as a list is
so that the monomials can appear in the same order as they are stored.
Using Plus
would revert to Mathematica’s default ordering.
For example
p = NCToNCPoly[1 + x**y**x - 2 x**x + z, vars]
NCPolyDisplay[p, vars]
returns
{x.y.x, z, -2 x.x, 1}
whereas
NCPolyToNC[p, vars]
would return
1 - 2 x^2 + z + x**y**x
in which the sorting of the monomials has been destroyed by
Plus
.
The monomials appear sorted in decreasing order from left to right,
with z
being the leading term in the above
example.
With NCPoly
the Mathematica command Sort
is
modified to sort lists of polynomials. For example
polys = NCToNCPoly[{x**x**x, 2 y**x - z, z, y**x - x**x}, vars];
ColumnForm[NCPolyDisplay[Sort[polys], vars]]
returns
{x.x.x}
{z}
{y.x, -x.x}
{2 y.x, -z}
Sort
produces a list of polynomials sorted in
ascending order based on their leading terms.
A larger class of polynomials in noncommutative variables is that of polynomials with noncommutative coefficients. Think of a polynomial with commutative coefficients in which certain variables are considered to be unknown, i.e. variables, where others are considered to be known, i.e. coefficients. For example, in many problems in systems and control the following expression
\[p(x) = a x + x a^T - x b x + c\]
is often seen as a polynomial in the noncommutative unknown
x
with known noncommutative coefficients a
,
b
, and c
. A typical problem is the
determination of a solution to the equation \(p(x) = 0\) or the inequality \(p(x) \succeq 0\).
The package NCPolynomial
handles such
polynomials with noncommutative coefficients. As with NCPoly
, the package provides the
commands NCToNCPolynomial
and NCPolynomialToNC
to
convert nc expressions back and forth between NCAlgebra
and
NCPolynomial
. For example
vars = {x}
p = NCToNCPolynomial[a**x + x**tp[a] - x**b**x + c, vars]
converts the polynomial a**x + x**tp[a] - x**b**x + c
from the standard NCAlgebra
format into an
NCPolynomial
object. The result in this case is the
NCPolynomial
object
NCPolynomial[c, <|{x} -> {{1, a, 1}, {1, 1, tp[a]}}, {x, x} -> {{-1, 1, b, 1}}|>, {x}]
Conversely the command NCPolynomialToNC
converts an
NCPolynomial
back into NCAlgebra
format. For
example
NCPolynomialToNC[p]
returns
c + a**x + x**tp[a] - x**b**x
An NCPolynomial
does store information about the
polynomial symbols and a list of variables is required only at the time
of creation of the NCPolynomial
object.
As with NCPoly
, operations on NCPolynomial
objects result on another NCPolynomial
object that is
always expanded. For example:
vars = {x,y}
1 + NCToNCPolynomial[x**y, vars] - 2 NCToNCPolynomial[y**x, vars]
returns
NCPolynomial[1, <|{y**x} -> {{-2, 1, 1}}, {x**y} -> {{1, 1, 1}}|>, {x, y}]
and
(1 + NCToNCPolynomial[x, vars]**NCToNCPolynomial[y, vars])^2
returns
NCPolynomial[1, <|{x**y**x**y} -> {{1, 1, 1}}, {x**y} -> {{2, 1, 1}}|>, {x, y}]
To see how much more efficient NCPolynomial
is when
compared with standard NCAlgebra
objects try
Table[Timing[(NCToNCPolynomial[x, vars])^i][[1]], {i, 0, 15, 5}]
would return
{0.000493, 0.003345, 0.005974, 0.013479, 0.018575}
As you can see, NCPolynomial
s are not as efficient as
NCPoly
s but still much more efficient than
NCAlgebra
polynomials.
NCPolynomials
do not support orderings but we
do provide the NCPSort
command that produces a list of
terms sorted by degree. For example
NCPSort[p]
returns
{c, a**x, x**tp[a], -x**b**x}
A useful feature of NCPolynomial
is the capability of
handling polynomial matrices. For example
mat1 = {{a**x + x**tp[a] + c**y + tp[y]**tp[c] - x**q**x, b**x},
{x**tp[b], 1}};
p1 = NCToNCPolynomial[mat1, {x, y}];
mat2 = {{1, x**tp[c]}, {c**x, 1}};
p2 = NCToNCPolynomial[mat2, {x, y}];
constructs NCPolynomial
objects representing the
polynomial matrices mat1
and mat2
. Verify
that
NCPolynomialToNC[p1**p2] - NCDot[mat1, mat2] // NCExpand
is zero as expected. Internally NCPolynomial
represents
a polynomial matrix by constructing matrix factors. For example the
representation of the matrix mat1
correspond to the factors
\[
\begin{aligned}
\begin{bmatrix}
a x + x a^T + c y + y^T c^T - x q x & b x \\
x b^T & 1
\end{bmatrix}
&=
\begin{bmatrix} 0 & 0 \\ 0 & 1 \end{bmatrix}
+
\begin{bmatrix} a \\ 0 \end{bmatrix}
x
\begin{bmatrix} 1 & 0 \end{bmatrix}
+
\begin{bmatrix} 1 \\ 0 \end{bmatrix}
x
\begin{bmatrix} a^T & 0 \end{bmatrix}
+
\begin{bmatrix} -1 \\ 0 \end{bmatrix}
x q x
\begin{bmatrix} 1 & 0 \end{bmatrix}
+ \\ & \qquad \quad
\begin{bmatrix} b \\ 0 \end{bmatrix}
x
\begin{bmatrix} 0 & 1 \end{bmatrix}
+
\begin{bmatrix} 0 \\ 1 \end{bmatrix}
x
\begin{bmatrix} b^T & 0 \end{bmatrix}
+
\begin{bmatrix} c \\ 0 \end{bmatrix}
y
\begin{bmatrix} 1 & 0 \end{bmatrix}
+
\begin{bmatrix} 1 \\ 0 \end{bmatrix}
y^T
\begin{bmatrix} c^T & 0 \end{bmatrix}
\end{aligned}
\]
See section linear polynomials for more features on linear polynomial matrices.
Another interesting class of nc polynomials is that of linear polynomials, which can be factored in the form:
\[s(x) = l (F \otimes x) r\]
where \(l\) and \(r\) are vectors with symbolic expressions and \(F\) is a numeric matrix. This functionality is in the package
<< NCSylvester`
Use the command NCToNCSylvester
to factor a
linear nc polynomial into the the above form. For example:
vars = {x, y};
expr = 1 + a**x + x**tp[a] - x + b**y**d + tp[d]**tp[y]**tp[b];
{const, lin} = NCToNCSylvester[expr, vars];
which returns
const = 1
and an Association
lin
containing the
factorization. For example
lin[x]
returns a list with the left and right vectors l
and
r
and the coefficient array F
.
{{1, a}, {1, a^T}, SparseArray[< 2 >, {2, 2}]}
which in this case is the matrix:
\[ \begin{bmatrix} -1 & 1\\ 1 & 0 \end{bmatrix} \]
and
lin[tp[y]]
returns
{{d^T}, {b^T}, SparseArray[< 1 >, {1, 1}]}
Note that transposes and adjoints are treated as independent variables.
Perhaps the most useful consequence of the above factorization is the
possibility of producing a linear polynomial which has the smallest
possible number of terms, as explaining in detail in (Oliveira 2012).
This is done automatically by NCSylvesterToNC
. For
example
vars = {x, y};
expr = a**x**c - a**x**d - a**y**c + a**y**d + b**x**c - b**x**d - b**y**c + b**y**d;
{const, lin} = NCToNCSylvester[expr, vars];
NCSylvesterToNC[{const, lin}]
produces:
(a + b)**x**(c - d) + (a + b)**y**(-c + d)
This factorization even works with linear matrix polynomials, and is used by the our semidefinite programming algorithm (see Chapter Semidefinite Programming) to factor linear matrix inequalities in the least possible number of terms. For example:
vars = {x};
expr = {{a**x + x**tp[a], b**x, tp[c]},
{x**tp[b], -1, tp[d]},
{c, d, -1}};
{const, lin} = NCToNCSylvester[expr, vars]
result in:
const = SparseArray[< 6 >, {3, 3}]
lin = <|x -> {{1, a, b}, {1, tp[a], tp[b]}, SparseArray[< 4 >, {9, 9}]}|>
See (Oliveira 2012) for details on the structure of the constant array \(F\) in this case.
Gröbner Basis are useful in the study of algebraic relations. The
package NCGBX
provides an implementation of a
noncommutative Gröbner Basis algorithm.
Starting with Version 6, the old
C++
version of our Groebner Basis Algorithm is no longer included.
If you want a living version of this chapter just run the notebook
NC/DEMOS/3_NCGroebnerBasis.nb
.
In order to load NCGBX
one types:
<< NC`
<< NCAlgebra`
<< NCGBX`
or simply
<< NCGBX`
if NC
and NCAlgebra
have already been
loaded.
Most commutative algebra packages contain commands based on Gröbner
Basis and uses of Gröbner Basis. For example, in Mathematica, the
Solve
command puts collections of equations in a
canonical form which, for simple collections, readily yields a
solution. Likewise, the Mathematica Eliminate
command tries
to convert a collection of \(m\)
polynomial equations (often called relations)
\[\begin{aligned} p_1(x_1,\ldots,x_n) &= 0 \\ p_2(x_1,\ldots,x_n) &= 0 \\ \vdots \quad & \quad \, \, \vdots \\ p_m(x_1,\ldots,x_n) &= 0 \end{aligned}\]
in variables \(x_1,x_2, \ldots x_n\) to a triangular form, that is a new collection of equations like
\[\begin{aligned} q_1(x_1) &= 0 \\ q_2(x_1,x_2) &= 0 \\ q_3(x_1,x_2) &= 0 \\ q_4(x_1,x_2,x_3)&=0 \\ \vdots \quad & \quad \, \, \vdots \\ q_{r}(x_1,\ldots,x_n) &= 0. \end{aligned}\]
Here the polynomials \(\{q_j: 1\le j\le k_2\}\) generate the same ideal that the polynomials \(\{p_j : 1\le j \le k_1\}\) generate. Therefore, the set of solutions to the collection of polynomial equations \(\{p_j=0: 1\le j\le k_1\}\) equals the set of solutions to the collection of polynomial equations \(\{q_j=0: 1\le j\le k_2\}\). This canonical form greatly simplifies the task of solving collections of polynomial equations by facilitating backsolving for \(x_j\) in terms of \(x_1,\ldots,x_{j-1}\).
Readers who would like to know more about Gröbner Basis may want to
read [CLS]. The noncommutatative version of the algorithm implemented by
NCGB
is loosely based on [Mora].
Before calculating a Gröbner Basis, one must declare which variables
will be used during the computation and must declare a monomial
order which can be done using SetMonomialOrder
as
in:
SetMonomialOrder[{a, b, c}, x];
The monomial ordering imposes a relationship between the variables which are used to sort the monomials in a polynomial. The ordering implied by the above command can be visualized using:
PrintMonomialOrder[]
which in this case prints:
\[a < b < c \ll x.\]
A user does not need to know theoretical background related to monomials orders. Indeed, as we shall see soon, in many engineering problems, it suffices to know which variables correspond to quantities which are known and which variables correspond to quantities which are unknown. If one is solving for a variable or desires to prove that a certain quantity is zero, then one would want to view that variable as unknown. In the above example, the symbol ‘\(\ll\)’ separate the knowns, \(a, b, c\), from the unknown, \(x\). For more details on orderings see Section Orderings.
Our goal is to calculate the Gröbner basis associated with the following relations (i.e. a list of polynomials):
\[\begin{aligned} a \, x \, a &= c, & a \, b &= 1, & b \, a &= 1. \end{aligned}\]
We shall use the word relation to mean a polynomial in noncommuting indeterminates. For example, if an analyst saw the equation \(A B = 1\) for matrices \(A\) and \(B\), then he might say that \(A\) and \(B\) satisfy the polynomial equation \(a\, b - 1 = 0\). An algebraist would say that \(a\, b - 1\) is a relation.
To calculate a Gröbner basis one defines a list of relations:
rels = {a**x**a - c, a**b - 1, b**a - 1}
and issues the command:
gb = NCMakeGB[rels, 10]
which should produce an output similar to:
* * * * * * * * * * * * * * * *
* * * NCPolyGroebner * * *
* * * * * * * * * * * * * * * *
* Monomial order: a < b < c << x
* Reduce and normalize initial set
> Initial set could not be reduced
* Computing initial set of obstructions
> MAJOR Iteration 1, 4 polys in the basis, 2 obstructions
> MAJOR Iteration 2, 5 polys in the basis, 2 obstructions
* Cleaning up...
* Found Groebner basis with 3 polynomials
* * * * * * * * * * * * * * * *
The number 10
in the call to NCMakeGB
is
very important because a finite GB may not exist. It instructs
NCMakeGB
to abort after 10
iterations if a GB
has not been found at that point.
The result of the above calculation is the list of relations in the form of a list of rules:
{x -> b**c**b, a**b -> 1, b**a -> 1}
Version 5: For efficiency,
NCMakeGB
returns a list of rules instead of a list of polynomials. The left-hand side of the rule is the leading monomial in the current order. This is incompatible with early versions, which returned a list of polynomials. You can recover the old behavior setting the optionReturnRules -> False
. This can be done in theNCMakeGB
command or globally throughSetOptions[ReturnRules -> False]
.
Our favorite format for displaying lists of relations is
ColumnForm
.
ColumnForm[gb]
which results in
x -> b**c**b
a**b -> 1
b**a -> 1
The rules in the output represent the relations in the GB
with the left-hand side of the rule being the leading monomial.
Replacing Rule
by Subtract
recovers the
relations but one would then loose the leading monomial as Mathematica
alphabetizes the resulting sum.
Someone not familiar with GB’s might find it instructive to note this output GB effectively solves the input equation
\[a \, x \, a - c = 0\]
under the assumptions that
\[\begin{aligned} b \, a - 1 &= 0, & a \, b - 1 & =0, \end{aligned}\]
that is \(a = b^{-1}\) and produces the expected result in the form of the relation:
\[ x = b \, c \, b. \]
For a slightly more challenging example consider the same monomial order as before:
SetMonomialOrder[{a, b, c}, x];
that is
\(a < b < c \ll x\)
and the relations:
\[\begin{aligned} a \, x - c &= 0, \\ a \, b \, a - a &= 0, \\ b \, a \, b - b &= 0, \end{aligned}\]
from which one can recognize the problem of solving the linear equation \(a \, x = c\) in terms of the pseudo-inverse \(b = a^\dagger\). The calculation:
gb = NCMakeGB[{a**x - c, a**b**a - a, b**a**b - b}, 10];
ColumnForm[gb]
finds the Gröbner basis:
a**x -> c
a**b**c -> c
a**b**a -> a
b**a**b -> b
In this case the Gröbner basis cannot quite solve the equations but it remarkably produces the necessary condition for existence of solutions:
\[ 0 = a \, b \, c - c = a \, a^\dagger c - c \]
that can be interpreted as \(c\) being in the range-space of \(a\).
Our goal now is to verify if it is possible to simplify the following expression:
expr = b^2**a^2 - a^2**b^2 + a**b**a
knowing that
\[ a \, b \, a = b \]
using Gröbner basis. With that in mind we set the ordering:
SetMonomialOrder[a,b];
and calculate the GB associated with the constraint:
rels = {a**b**a - b};
rules = NCMakeGB[rels, 10];
ColumnForm[rules]
which produces the output
* * * * * * * * * * * * * * * *
* * * NCPolyGroebner * * *
* * * * * * * * * * * * * * * *
* Monomial order: a << b
* Reduce and normalize initial set
> Initial set could not be reduced
* Computing initial set of obstructions
> MAJOR Iteration 1, 2 polys in the basis, 1 obstructions
* Cleaning up...
* Found Groebner basis with 2 polynomials
* * * * * * * * * * * * * * * *
and the associated GB
a**b**a -> b
b^2**a -> a**b^2
The GB revealed another relationship that must hold true if \(a \, b \, a = b\). One can use these
relationships to simplify the original expression using
NCReplaceRepeated
as in
NCReplaceRepeated[expr, rules]
which simplifies expr
into b
.
An alternative, since we are working with polynomials, is to use NCReduce
as in
vars = GetMonomialOrder[];
NCReduce[expr, NCRuleToPoly[rules], vars]
Note that NCReduce
needs a list of variables to be used
as the desired ordering, which we obtain from the current ordering using
GetMonomialOrder
, and that
the rules in rules
need to be converted to polynomials,
which we do using NCRuleToPoly
.
Having seen how polynomial relations can be interpreted as rules, consider the expression
expr = b^2**a^2 + a^2**b^3 + a**b**a
and the polynomial relations
rels = {a**b**a - b , b^2 - a + b}
With respect to the monomial ordering
SetMonomialOrder[a, b];
we can interpret these relations as the rules
vars = GetMonomialOrder[];
(rules = NCToRule[rels, vars]) // ColumnForm
that is
a**b**a -> b
b^2 -> a - b
Note how we used GetMonomialOrder
to obtain the
list of variables vars
corresponding to the ordering
\[ a \ll b\]
and NCToRule
to convert the
polynomial relations into rules.
We can then apply these rules by using one of the
NCReplace
functions, for example
NCExpandReplaceRepeated[expr, rules]
which produces
b + a^2**b + a^3**b - b**a^2
Note that we have made use of the new function NCExpandReplaceRepeated
,
to automate the tedious cycles of expansion and substitution.
Alternatively, one can use NCReduce
to perform the same
substitution. NCReduce
takes in polynomial, instead of
rules, and a list of variables. For example,
NCReduce[expr, rels, vars]
produces
a^2**b + a^3**b - b**a^2 + a**b**a
which is the result of applying the rules only to the leading
monomial of expr
. If you want substitutions to be applied
to all monomials then set the option Complete
to
True
, as in
NCReduce[expr, rels, vars, Complete -> True]
This produces
b + a^2**b + a^3**b - b**a^2
which is the same result as above.
But, of course, by now one should wonder whether the above polynomial relations could imply other simplifications, which we seek to find out by running our Gröbner basis algorithm:
rules = NCMakeGB[rels, 4];
ColumnForm[rules]
that discovers the following additional relations
b^2 -> a - b
b**a -> a**b
a^2**b -> b
a^3 -> a
When used for simplification,
NCExpandReplaceRepeated[expr, rules]
reduces the original expression to the even simpler form
b + a**b
As before, rule substitution could also be performed by NCReduce
as in
NCReduce[expr, NCRuleToPoly[rules], vars]
that, in this case, leads to the same result as above without the
need to recourse to the Complete
flag.
The algorithm implemented by NCGB
always produces a
Gröbner Basis with the minimal possible number of polynomials
at a given iteration. However, such polynomials are not necessarily the
“simplest” possible polynomials; called the reduced Gröbner
Basis. The reduced Gröbner Basis is unique given the relations
and the monomial ordering. Consider for example the following monomial
ordering
SetMonomialOrder[x, y]
and the relations
rels = {x^3 - 2 x**y, x^2**y - 2 y^2 + x}
for which
NCMakeGB[rels, ReduceBasis -> False] // ColumnForm
produces the minimal Gröbner Basis
x^2->0
x**y->x^3/2
y**x->x**y
y^2->x/2+x^2**y/2
but
NCMakeGB[rels, ReduceBasis -> True] // ColumnForm
returns the reduced Gröbner Basis
x^2->0
x**y->0
y**x->0
y^2->x/2
in which not only the leading monomials but also all lower-order
monomials have been reduced by the basis’ leading monomials. The option
ReduceBasis
is set to True
by default.
It is often desirable to simplify expressions involving inverses of noncommutative expressions. One challenge is to recognize identities implied by the existence of certain inverses. For example, that the expression
expr = x**inv[1 - x] - inv[1 - x]**x
is equivalent to \(0\). One can use a nc Gröbner basis for that task. Consider for instance the ordering
\[ x \ll (1-x)^{-1} \]
implied by the command:
SetMonomialOrder[x, inv[1-x]]
This ordering encodes the following precise idea of what we mean by simple versus complicated: it formally corresponds to specifying that \(x\) is simpler than \((1-x)^{-1}\), which might sits well with one’s intuition.
Now consider the following command:
(rules = NCMakeGB[{}, 3]) // ColumnForm
which produces the output
* * * * * * * * * * * * * * * *
* * * NCPolyGroebner * * *
* * * * * * * * * * * * * * * *
* Monomial order: x << inv[x] << inv[1 - x]
* Reduce and normalize initial set
> Initial set could not be reduced
* Computing initial set of obstructions
> MAJOR Iteration 1, 6 polys in the basis, 6 obstructions
* Cleaning up...
* Found Groebner basis with 6 polynomials
* * * * * * * * * * * * * * * *
and results in the rules:
x**inv[1 - x] -> -1 + inv[1 - x],
inv[1-x]**x -> -1 + inv[1-x],
As in the previous example, the GB revealed new relationships that
must hold true if \(1- x\) is
invertible, and one can use this relationship to simplify the
original expression using NCReplaceRepeated
as in:
NCReplaceRepeated[expr, rules]
The above command results in 0
, as one would hope.
For a more challenging example consider the identity:
\[\left (1 - x - y (1 - x)^{-1} y \right )^{-1} = \frac{1}{2} (1 - x - y)^{-1} + \frac{1}{2} (1 - x + y)^{-1}\]
One can verify that the rule based command NCSimplifyRational fails to simplify the expression:
expr = inv[1 - x - y**inv[1 - x]**y] - 1/2 (inv[1 - x + y] + inv[1 - x - y])
NCSimplifyRational[expr]
We set the monomial ordering and calculate the Gröbner basis
SetMonomialOrder[x, y, inv[1-x], inv[1-x+y], inv[1-x-y], inv[1-x-y**inv[1-x]**y]];
(rules = NCMakeGB[{}, 3]) // ColumnForm
based on the rational involved in the original expression. The result is the nc GB:
inv[1-x-y**inv[1-x]**y] -> (1/2)inv[1-x-y]+(1/2)inv[1-x+y]
x**inv[1-x] -> -1+inv[1-x]
y**inv[1-x+y] -> 1-inv[1-x+y]+x**inv[1-x+y]
y**inv[1-x-y] -> -1+inv[1-x-y]-x**inv[1-x-y]
inv[1-x]**x -> -1+inv[1-x]
inv[1-x+y]**y -> 1-inv[1-x+y]+inv[1-x+y]**x
inv[1-x-y]**y -> -1+inv[1-x-y]-inv[1-x-y]**x
inv[1-x+y]**x**inv[1-x-y] -> -(1/2)inv[1-x-y]-(1/2)inv[1-x+y]+inv[1-x+y]**inv[1-x-y]
inv[1-x-y]**x**inv[1-x+y] -> -(1/2)inv[1-x-y]-(1/2)inv[1-x+y]+inv[1-x-y]**inv[1-x+y]
which successfully simplifies the original expression using:
NCExpandReplaceRepeated[expr, rules]
resulting in 0
.
See also Advanced Processing of Rational Expressions for more details on the lower level handling of rational expressions.
The simplification process described above is automated in the function NCGBSimplifyRational.
For example, calls to
expr = x**inv[1 - x] - inv[1 - x]**x
NCGBSimplifyRational[expr]
or
expr = inv[1 - x - y**inv[1 - x]**y] - 1/2 (inv[1 - x + y] + inv[1 - x - y])
NCGBSimplifyRational[expr]
both result in 0
.
As seen above, one needs to declare a monomial order before making a Gröbner Basis. There are various monomial orders which can be used when computing Gröbner Basis. The most common are lexicographic and graded lexicographic orders. We consider also multi-graded lexicographic orders.
Lexicographic and multi-graded lexicographic orders are examples of elimination orderings. An elimination ordering is an ordering which is used for solving for some of the variables in terms of others.
We now discuss each of these types of orders.
To impose lexicographic order, say \(a\ll b\ll x\ll y\) on \(a\), \(b\), \(x\) and \(y\), one types
SetMonomialOrder[a,b,x,y];
or, equivalently
SetMonomialOrder[{a},{b},{x},{y}];
This order is useful for attempting to solve for \(y\) in terms of \(a\), \(b\) and \(x\), since the highest priority of the GB algorithm is to produce polynomials which do not contain \(y\). If producing high order polynomials is a consequence of this fanaticism so be it. Unlike graded orders, lex orders pay little attention to the degree of terms. Likewise its second highest priority is to eliminate \(x\).
Once this order is set, one can use all of the commands in the preceding section in exactly the same form.
We now give a simple example how one can solve for \(y\) given that \(a\),\(b\),\(x\) and \(y\) satisfy the equations:
\[\begin{aligned} -b\, x + x\, y \, a + x\, b \, a \, a &= 0 \\ x \, a-1&=0 \\ a\, x-1&=0 \end{aligned}\]
The command
NCMakeGB[{-b**x+x**y**a+x**b**a**a, x**a-1, a**x-1},4] // ColumnForm
produces the Gröbner basis:
y -> -b**a + a**b**x^2
a**x -> 1
x**a -> 1
after one iteration.
Now, we change the ordering to
SetMonomialOrder[y,x,b,a];
and run the same NCMakeGB
as above:
NCMakeGB[{-b**x+x**y**a+x**b**a**a, x**a-1, a**x-1},4] // ColumnForm
which, this time, results in
b**x^3 -> x**b+x**y**x
x**a -> 1
a**x -> 1
x**b**a -> b**x^2 - x**y
a**b**x^2 -> y+b**a
x**b^2**a -> -x**b**y+b**x^2**b**x^2 -
x**y**b**x^2
b**a^2 -> -y**a+a**b**x
b**a**b**a -> -y^2 - b**a**y - y**b**a+
a**b**x**b**x^2
b**a**b^2**a -> -y**b**y - y**b^2**a -
y^2**b**x^2 - b**a**b**y - b**a**y**b**x^2+
a**b**x**b**x^2**b**x^2
which is not a Gröbner basis since the algorithm was interrupted at 4 iterations. Note the presence of the rule
a**b**x**x -> y+b**a
which shows that the ordering is not set up to solve for \(y\) in terms of the other variables in the sense that \(y\) is not on the left hand side of this rule (but a human could easily solve for \(y\) using this rule). Also the algorithm created a number of other relations which involved \(y\).
To impose graded lexicographic order, say \(a< b< x< y\) on \(a\), \(b\), \(x\) and \(y\), one types
SetMonomialOrder[{a,b,x,y}];
This ordering puts high degree monomials high in the ordering. Thus it tries to decrease the total degree of expressions. A call to
NCMakeGB[{-b**x+x**y**a+x**b**a**a, x**a-1, a**x-1},4,ReduceBasis->True] // ColumnForm
now produces
a**x -> 1,
x**a -> 1,
b**a^2 -> -y**a+a**b**x,
x**b**a -> b**x^2 - x**y,
a**b**x^2 -> y+b**a,
b**x^3 -> x**b+x**y**x,
a**b**x**b**x^2 -> y^2+b**a**y+y**b**a+b**a**b**a,
b**x^2**b**x^2 -> x**b**y+x**b^2**a+x**y**b**x^2,
a**b**x**b**x**b**x^2 -> y^3+b**a**y^2+y^2**b**a+y**b**a**y+
b**a**b**a**y+b**a**y**b**a+
y**b**a**b**a+b**a**b**a**b**a,
b**x^2**b**x**b**x^2 -> x**b**y^2+x**b^2**a**y+x**b**y**b**a+
x**b^2**a**b**a+x**y**b**x**b**x^2
which again fails to be a Gröbner basis and does not eliminate \(y\). Instead, it tries to decrease the total degree of expressions involving \(a\), \(b\), \(x\), and \(y\).
There are other useful monomial orders which one can use other than graded lex and lex. Another type of order is what we call multigraded lex and is a mixture of graded lex and lex order. To impose multi-graded lexicographic order, say \(a< b< x\ll y\) on \(a\), \(b\), \(x\) and \(y\), one types
SetMonomialOrder[{a,b,x},y];
which separates \(y\) from the remaining variables. This time, a call to
NCMakeGB[{-b**x+x**y**a+x**b**a**a, x**a-1, a**x-1},4,ReduceBasis->True] // ColumnForm
yields once again
y -> -b**a+a**b**x^2
a**x -> 1
x**a -> 1
which not only eliminates \(y\) but is also Gröbner basis, calculated after one iteration.
For an intuitive idea of why multigraded lex is helpful, we think of \(a\), \(b\), and \(x\) as corresponding to variables in some engineering problem which represent quantities which are known and \(y\) to be unknown. The fact that \(a\), \(b\) and \(x\) are in the top level indicates that we are very interested in solving for \(y\) in terms of \(a\), \(b\), and \(x\), but are not willing to solve for, say \(x\), in terms of expressions involving \(y\).
This situation is so common that we provide the commands
SetKnowns
and SetUnknowns
. The above ordering
would be obtained after setting
SetKnowns[a,b,x];
SetUnknowns[y];
This is a type of problem known as a matrix completion problem. This particular one was suggested by Hugo Woerdeman. We are grateful to him for discussions.
Problem: Given matrices \(a\), \(b\), \(c\), and \(d\), we wish to determine under what conditions there exists matrices x, y, z, and w such that the block matrices
\[\begin{bmatrix} a & x \\ y & b \end{bmatrix} \qquad \begin{bmatrix} w & c \\ d & z \end{bmatrix}\]
are inverses of each other. Also, we wish to find formulas for \(x\), \(y\), \(z\), and \(w\).
This problem was solved in a paper by W.W. Barrett, C.R. Johnson, M. E. Lundquist and H. Woerderman [BJLW] where they showed it splits into several cases depending upon which of \(a\), \(b\), \(c\) and \(d\) are invertible. In our example, we assume that \(a\), \(b\), \(c\) and \(d\) are invertible and discover the result which they obtain in this case.
First we set the matrices \(a\), \(b\), \(c\), and \(d\) and their inverses as knowns and \(x\), \(y\), \(w\), and \(z\) as unknowns:
SetKnowns[a, inv[a], b, inv[b], c, inv[c], d, inv[d]];
SetUnknowns[{z}, {x, y, w}];
Note that the graded ordering of the unknowns means that we care more about solving for \(x\), \(y\) and \(w\) than for \(z\).
Then we define the relations we are interested in, which are obtained after multiplying the two block matrices on both sides and equating to identity
A = {{a, x}, {y, b}}
B = {{w, c}, {d, z}}
rels = {
NCDot[A, B] - IdentityMatrix[2],
NCDot[B, A] - IdentityMatrix[2]
} // Flatten
We use Flatten
to reduce the matrix relations to a
simple list of relations. The resulting relations in this case are:
rels = {-1+a**w+x**d, a**c+x**z, b**d+y**w, -1+b**z+y**c,
-1+c**y+w**a, c**b+w**x, d**a+z**y, -1+d**x+z**b}
After running
NCMakeGB[rels, 8] // ColumnForm
we obtain the Gröbner basis:
x -> inv[d]-inv[d]**z**b
y -> inv[c]-b**z**inv[c]
w -> inv[a]**inv[d]**z**b**d
z**b**z -> z+d**a**c
c**b**z**inv[c]**inv[a] -> inv[a]**inv[d]**z**b**d
inv[c]**inv[a]**inv[d]**z**b -> b**z**inv[c]**inv[a]**inv[d]
inv[d]**z**b**d**a -> a**c**b**z**inv[c]
z**b**d**a**c -> d**a**c**b**z
z**inv[c]**inv[a]**inv[d]**inv[b] -> inv[b]**inv[c]**inv[a]**inv[d]**z
z**inv[c]**inv[a]**inv[d]**z -> inv[b]+inv[b]**inv[c]**inv[a]**inv[d]**z
d**a**c**b**z**inv[c] -> z**b**d**a
after seven iterations. The first four relations
\[\begin{aligned} x &= d^{-1}-d^{-1} \, z \, b \\ y &= c^{-1}-b \, z \, c^{-1} \\ w &= a^{-1} \, d^{-1} \, z \, b \, d \\ z \, b \, z &= z + d \, a \, c \end{aligned}\]
are the solutions we are looking for, which states that one can find \(x\), \(y\), \(z\), and \(w\) such that the matrices above are inverses of each other if and only if \(z \, b \, z = z + d \, a \, c\). The first three relations gives formulas for \(x\), \(y\) and \(w\) in terms of \(z\).
A variety of scenarios can be quickly investigated under different assumptions. For example, say that \(c\) is not invertible. Is it still possible to solve the problem? One solution is obtained with the ordering implied by
SetKnowns[a, inv[a], b, inv[b], c, d, inv[d]];
SetUnknowns[{y}, {z, w, x}];
In this case
NCMakeGB[rels, 8] // ColumnForm
produces the Gröbner basis:
z -> inv[b]-inv[b]**y**c
w -> inv[a]-c**y**inv[a]
x -> a**c**y**inv[a]**inv[d]
y**c**y -> y+b**d**a
c**y**inv[a]**inv[d]**inv[b] -> inv[a]**inv[d]**inv[b]**y**c
d**a**c**y**inv[a] -> inv[b]**y**c**b**d
inv[d]**inv[b]**y**c**b -> a**c**y**inv[a]**inv[d]
y**c**b**d**a -> b**d**a**c**y
y**inv[a]**inv[d]**inv[b]**y**c -> 1+y**inv[a]**inv[d]**inv[b]
after five iterations. Once again, the first four relations
\[\begin{aligned} z &= b^{-1}-b^{-1} \, y \, c \\ w &= a^{-1}-c \, y \, a^{-1} \\ x &= a \, c \, y \, a^{-1} \, d^{-1} \\ y \, c \, y &= y+b \, d \, a \end{aligned}\]
provide formulas, this time for \(z\), \(w\), and \(z\) in terms of \(y\) satisfying \(y \, c \, y = y+b \, d \, a\). Note that these formulas do not involve \(c^{-1}\) since \(c\) is no longer assumed invertible.
Consider once again the task of simplifying the nc rational expression
expr = inv[1-x-y**inv[1-x]**y] - 1/2 (inv[1-x+y]+inv[1-x-y])
considered before in Simplifying Rational Expressions. We will use this expression to illustrate how nc rational expressions can be manipulated at a lower level, giving advanced users more control of the conversion process to and from the internal NCPoly representation.
The key functionality is provided by the functions NCMonomialOrder and NCRationalToNCPoly. NCMonomialOrder provides the same functionality as SetMonomialOrder, but instead of setting the ordering globally, it returns an array representing the ordering. For example,
order = NCMonomialOrder[x, y]
produces the array
{{x},{y}}
which represents the ordering \(x \ll y\).
With a (preliminary) ordering in hand one can invoke NCRationalToNCPoly as in
{rels, vars, rules, labels} = NCRationalToNCPoly[expr, order];
This function produces four lists as outputs:
The first, rels
, is a list of NCPoly
objects representing the original expression and some additional
relations that were automatically generated. We will inspect
rels
later.
The second, vars
, is a list of variables that represent
the ordering used in the construction of the polynomials in
rel
. In this example,
vars = {{x}, {y}, {rat135, rat136, rat137, rat138}}
which corresponds to the ordering
\[x \ll y \ll rat135 < rat136 < rat137 < rat138\]
In this case, the variables rat135
, rat136
,
rat137
, rat138
were automatically created and
assigned an ordering by NCRationalToNCPoly
.
It is unlikely that your variables will be assigned the same suffixes
135
through138
. Those suffixes were chosen by Mathematica to make these variables unique in the global context, and hence depend on a multitude of factors that can only be determined at the time of execution.
As with SetMonomialOrder and NCMakeGB, rational terms can be explicitly added to
the ordering for finer control. For example, calling
NCRationalToNCPoly
with
order = NCMonomialOrder[x,y,inv[1-x],{inv[1-x+y],inv[1-x-y]}];
would have produced
vars = {{x}, {y}, {rat135}, {rat136, rat137}, {rat138}}
which corresponds to the ordering
\[x \ll y \ll rat135 \ll rat136 < rat137 \ll rat138\]
The third is a list of rules relating the new variables with rational
terms appearing in the original expression expr
. In this
example
rules = {rat135 -> inv[1-x], rat136 -> inv[1-x-y],
rat137 -> inv[1-x+y], rat138 -> inv[1-x-y**rat135**y]}
Finally, the fourth is a list of labels that is used for printing or displaying
labels = {{x}, {y}, {inv[1-x], inv[1-x-y], inv[1-x+y],
inv[1-x-y**inv[1-x]**y]}}
as used by NCMakeGB
and NCPolyDisplay
.
The relations in rels
can be visualized by using
NCPolyToNC
. For example,
NCPolyToNC[#, vars] & /@ rels // ColumnForm
produces
rat138-rat136/2-rat137/2
rat135-rat135**x-1
rat135-x**rat135-1
rat136-rat136**x-rat136**y-1
rat136-x**rat136-y**rat136-1
rat137-rat137**x+rat137**y-1
rat137-x**rat137+y**rat137-1
rat138-rat138**x-rat138**y**rat135**y-1
rat138-x**rat138-y**rat135**y**rat138-1
The first entry is simply the original expr
in which
every rational expression has been substituted by a new variable. The
same could be obtained by applying reverse rules
and
applying it repeatedly
NCReplaceRepeated[expr, Reverse /@ rules]
The remaining entries are polynomials encoding the rational expressions that have been substituted by new variables. For example, the first two additional relations,
rat135-rat135**x-1
rat135-x**rat135-1
correspond to the assertion that rat153 == inv[1-x]
, and
so on.
Equipped with a set of polynomial relations encoding the rational
expression expr
one can calculated a Gröebner basis by
calling the low-level implementation NCPolyGroebner. In this example, calling
{basis, tree} = NCPolyGroebner[Rest[rels], 4, Labels -> labels];
produces the output
* * * * * * * * * * * * * * * *
* * * NCPolyGroebner * * *
* * * * * * * * * * * * * * * *
* Monomial order: x<<y<<inv[1-x]<inv[1-x-y]<inv[1-x+y]<inv[1-x-y**inv[1-x]**y]
* Reduce and normalize initial set
> Initial set could not be reduced
* Computing initial set of obstructions
> MAJOR Iteration 1, 10 polys in the basis, 12 obstructions
> MAJOR Iteration 2, 10 polys in the basis, 6 obstructions
> Found Groebner basis with 9 polynomials
* * * * * * * * * * * * * * * *
Note the use of labels
to pretty print the monomial
ordering. The resulting basis can be visualized once again using the
labels
NCPolyToNC[#, labels] & /@ basis // ColumnForm
which produces
1-inv[1-x]+inv[1-x]**x
1-inv[1-x]+x**inv[1-x]
1-inv[1-x-y]+inv[1-x-y]**x+inv[1-x-y]**y
1-inv[1-x-y]+x**inv[1-x-y]+y**inv[1-x-y]
-1+inv[1-x+y]-inv[1-x+y]**x+inv[1-x+y]**y
-1+inv[1-x+y]-x**inv[1-x+y]+y**inv[1-x+y]
1/2 inv[1-x-y]+1/2 inv[1-x+y]-inv[1-x+y]**inv[1-x-y]+inv[1-x+y]**x**inv[1-x-y]
1/2 inv[1-x-y]+1/2 inv[1-x+y]-inv[1-x-y]**inv[1-x+y]+inv[1-x-y]**x**inv[1-x+y]}
-1/2 inv[1-x-y]-1/2 inv[1-x+y]+inv[1-x-y**inv[1-x]**y]
The original expression expr
can then be
reduced by the above basis by calling
NCPolyReduce[rels[[1]], basis]
which produces 0
, as expected.
The above process is conveniently automated by NCMakeGB, but advanced users might still want to
take advantage of the increased speed of directly processing
NCPoly
s by manually performing the conversion from a
rational statement to a polynomial statement. See the detailed
documentation of the package NCPoly for more
details.
If you want a living version of this chapter just run the notebook
NC/DEMOS/4_SemidefiniteProgramming.nb
.
There are two different packages for solving semidefinite programs:
SDP
provides a template
algorithm that can be customized to solve semidefinite programs with
special structure. Users can provide their own functions to evaluate the
primal and dual constraints and the associated Newton system. A built in
solver along conventional lines, working on vector variables, is
provided by default. It does not require NCAlgebra to run.
NCSDP
coordinates with
NCAlgebra to handle matrix variables, allowing constraints, etc, to be
entered directly as noncommutative expressions.
The package NCSDP allows the symbolic manipulation and numeric solution of semidefinite programs.
After loading NCAlgebra, the package NCSDP must be loaded using:
<< NCSDP`
Semidefinite programs consist of symbolic noncommutative expressions representing inequalities and a list of rules for data replacement. For example the semidefinite program: \[ \begin{aligned} \min_Y \quad & <I,Y> \\ \text{s.t.} \quad & A Y + Y A^T + I \preceq 0 \\ & Y \succeq 0 \end{aligned} \] can be solved by defining the noncommutative expressions
SNC[a, y];
obj = {-1};
ineqs = {a ** y + y ** tp[a] + 1, -y};
The inequalities are stored in the list ineqs
in the
form of noncommutative linear polyonomials in the variable
y
and the objective function constains the symbolic
coefficients of the inner product, in this case -1
. The
reason for the negative signs in the objective as well as in the second
inequality is that semidefinite programs are expected to be cast in the
following canonical form: \[
\begin{aligned}
\max_y \quad & <b,y> \\
\text{s.t.} \quad & f(y) \preceq 0
\end{aligned}
\] or, equivalently: \[
\begin{aligned}
\max_y \quad & <b,y> \\
\text{s.t.} \quad & f(y) + s = 0, \quad s \succeq 0
\end{aligned}
\]
Semidefinite programs can be visualized using NCSDPForm
as in:
vars = {y};
NCSDPForm[ineqs, vars, obj]
The above commands produce a formatted output similar to the ones shown above.
In order to obtaining a numerical solution for an instance of the above semidefinite program one must provide a list of rules for data substitution. For example:
A = {{0, 1}, {-1, -2}};
data = {a -> A};
Equipped with the above list of rules representing a problem instance
one can load SDPSylvester
and use
NCSDP
to create a problem instance as follows:
{abc, rules} = NCSDP[ineqs, vars, obj, data];
The resulting abc
and rules
objects are
used for calculating the numerical solution using SDPSolve
. The command:
<< SDPSylvester`
{Y, X, S, flags} = SDPSolve[abc, rules];
produces an output like the folowing:
Problem data:
* Dimensions (total):
- Variables = 4
- Inequalities = 2
* Dimensions (detail):
- Variables = {{2,2}}
- Inequalities = {2,2}
Method:
* Method = PredictorCorrector
* Search direction = NT
Precision:
* Gap tolerance = 1.*10^(-9)
* Feasibility tolerance = 1.*10^(-6)
* Rationalize iterates = False
Other options:
* Debug level = 0
K <B, Y> mu theta/tau alpha |X S|2 |X S|oo |A* X-B| |A Y+S-C|
-------------------------------------------------------------------------------------------
1 1.638e+00 1.846e-01 2.371e-01 8.299e-01 1.135e+00 9.968e-01 9.868e-16 2.662e-16
2 1.950e+00 1.971e-02 2.014e-02 8.990e-01 1.512e+00 9.138e-01 2.218e-15 2.937e-16
3 1.995e+00 1.976e-03 1.980e-03 8.998e-01 1.487e+00 9.091e-01 1.926e-15 3.119e-16
4 2.000e+00 9.826e-07 9.826e-07 9.995e-01 1.485e+00 9.047e-01 8.581e-15 2.312e-16
5 2.000e+00 4.913e-10 4.913e-10 9.995e-01 1.485e+00 9.047e-01 1.174e-14 4.786e-16
-------------------------------------------------------------------------------------------
* Primal solution is not strictly feasible but is within tolerance
(0 <= max eig(A* Y - C) = 8.06666*10^-10 < 1.*10^-6 )
* Dual solution is within tolerance
(|| A X - B || = 1.96528*10^-9 < 1.*10^-6)
* Feasibility radius = 0.999998
(should be less than 1 when feasible)
The output variables Y
and S
are the
primal solutions and X
is the dual
solution.
A symbolic dual problem can be calculated easily using NCSDPDual
:
{dIneqs, dVars, dObj} = NCSDPDual[ineqs, vars, obj];
The dual program for the example problem above is: \[
\begin{aligned}
\max_x \quad & <c,x> \\
\text{s.t.} \quad & f^*(x) + b = 0, \quad x \succeq 0
\end{aligned}
\] In the case of the above problem the dual program is \[
\begin{aligned}
\max_{X_1, X_2} \quad & <I,X_1> \\
\text{s.t.} \quad & A^T X_1 + X_1 A -X_2 - I = 0 \\
& X_1 \succeq 0, \\
& X_2 \succeq 0
\end{aligned}
\] which can be visualized using NCSDPDualForm
using:
NCSDPDualForm[dIneqs, dVars, dObj]
The package SDP provides a crude and not
very efficient way to define and solve semidefinite programs in standard
form, that is vectorized. You do not need to load NCAlgebra
if you just want to use the semidefinite program solver. But you still
need to load NC
as in:
<< NC`
<< SDP`
Semidefinite programs are optimization problems of the form: \[ \begin{aligned} \max_{y, S} \quad & b^T y \\ \text{s.t.} \quad & A y + S = c \\ & S \succeq 0 \end{aligned} \] where \(S\) is a symmetric positive semidefinite matrix and \(y\) is a vector of decision variables.
A user can input the problem data, the triplet \((A, b, c)\), or use the following convenient methods for producing data in the proper format.
For example, problems can be stated as: \[ \begin{aligned} \min_y \quad & f(y), \\ \text{s.t.} \quad & G(y) \succeq 0 \end{aligned} \] where \(f(y)\) and \(G(y)\) are affine functions of the vector of variables \(y\).
Here is a simple example:
y = {y0, y1, y2};
f = y2;
G = {y0 - 2, {{y1, y0}, {y0, 1}}, {{y2, y1}, {y1, 1}}};
The list of constraints in G
is to be interpreted as:
\[
\begin{aligned}
y_0 - 2 \geq 0, \\
\begin{bmatrix} y_1 & y_0 \\ y_0 & 1 \end{bmatrix} \succeq 0,
\\
\begin{bmatrix} y_2 & y_1 \\ y_1 & 1 \end{bmatrix} \succeq 0.
\end{aligned}
\] The function SDPMatrices
convert the above
symbolic problem into numerical data that can be used to solve an
SDP.
abc = SDPMatrices[f, G, y]
All required data, that is \(A\),
\(b\), and \(c\), is stored in the variable
abc
as Mathematica’s sparse matrices. Their contents can be
revealed using the Mathematica command Normal
.
Normal[abc]
The resulting SDP is solved using SDPSolve
:
{Y, X, S, flags} = SDPSolve[abc];
The variables Y
and S
are the
primal solutions and X
is the dual
solution. Detailed information on the computed solution is found in the
variable flags
.
The package SDP
is built so as to be easily overloaded
with more efficient or more structure functions. See for example SDPFlat and SDPSylvester.
If you want a living version of this chapter just run the notebook
NC/DEMOS/5_PrettyOutput.nb
.
NCAlgebra
comes with several utilities for beautifying
expressions which are output. NCTeXForm
converts NC
expressions into . NCTeX
goes a
step further and compiles the results expression in and produces a PDF
that can be embedded in notebooks of used on its own.
In a Mathematica notebook session the package NCOutput can be used to control how nc
expressions are displayed. NCOutput
does not alter the
internal representation of nc expressions, just the way they are
displayed on the screen.
The function NCSetOutput can be used to set the display options. For example:
NCSetOutput[tp -> False, inv -> True];
makes the expression
expr = inv[tp[a] + b]
be displayed as \[ (\mathrm{tp[a]} + \mathrm{b})^{-1} \] Conversely
NCSetOutput[tp -> True, inv -> False];
makes expr
be displayed as \[
\mathrm{inv}\mathrm{[}\mathrm{a}^T + \mathrm{b}\mathrm{]}
\]
The default settings are
NCSetOutput[tp -> True, inv -> True];
which makes expr
be displayed as \[
(\mathrm{a}^\mathrm{T} + \mathrm{b})^{-1}
\] The complete set of options and their default values are:
NonCommutativeMultiply
(False
): If
True
x**y
is displayed as ‘\(\mathrm{x} \bullet \mathrm{y}\)’;tp
(True
): If True
tp[x]
is displayed as ‘\(\mathrm{x}^\mathrm{T}\)’;inv
(True
): If True
inv[x]
is displayed as ‘\(\mathrm{x}^{-1}\)’;aj
(True
): If True
aj[x]
is displayed as ‘\(\mathrm{x}^*\)’;co
(True
): If True
co[x]
is displayed as ‘\(\bar{\mathrm{x}}\)’;rt
(True
): If True
rt[x]
is displayed as ‘\(\mathrm{x}^{1/2}\)’.The special symbol All
can be used to set all options to
True
or False
, as in
NCSetOutput[All -> True];
You can load NCTeX using the following command
<< NC`
<< NCTeX`
NCTeX
does not need NCAlgebra
to work. You
may want to use it even when not using NCAlgebra. It uses NCRun
, which is a replacement for
Mathematica’s Run command to run pdflatex
,
latex
, divps
, etc.
WARNING: Mathematica does not come with LaTeX,
dvips, etc. The package NCTeX
does not install these
programs but rather assumes that they have been previously installed and
are available at the user’s standard shell. Use the Verbose
option to troubleshoot installation
problems.
With NCTeX
loaded you simply type
NCTeX[expr]
and your expression will be converted to a PDF
image which, by default, appears in your notebook after being processed
by LaTeX
. See options for
information on how to change this behavior to display the PDF on a
separate window.
For example:
expr = 1 + Sin[x + (y - z)/Sqrt[2]];
NCTeX[expr]
produces
\(1 + \sin \left ( x + \frac{y - z}{\sqrt{2}} \right )\)
If NCAlgebra
is not loaded then NCTeX
uses
the built in TeXForm
to produce the LaTeX expressions. If
NCAlgebra
is loaded, NCTeXForm
is used. See NCTeXForm for details.
Here is another example:
expr = {{1 + Sin[x + (y - z)/2 Sqrt[2]], x/y}, {z, n Sqrt[5]}};
NCTeX[expr]
that produces
\(\left( \begin{array}{cc} \sin \left(x+\frac{y-z}{\sqrt{2}}\right)+1 & \frac{x}{y} \\ z & \sqrt{5} n \\ \end{array} \right)\)
In some cases Mathematica will have difficulty displaying certain PDF
files. When this happens NCTeX
will span a PDF viewer so
that you can look at the formula. If your PDF viewer does not pop up
automatically you can force it by passing the following option to
NCTeX
:
expr = {{1 + Sin[x + (y - z)/2 Sqrt[2]], x/y}, {z, n Sqrt[5]}};
NCTeX[exp, DisplayPDF -> True]
Here is another example were the current version of Mathematica fails to import the PDF:
expr = Table[x^i y^(-j) , {i, 0, 10}, {j, 0, 30}];
NCTeX[expr, DisplayPDF -> True]
You can also suppress Mathematica from importing the PDF altogether as well. This and other options are covered in detail in the next section.
The following command:
expr = {{1 + Sin[x + (y - z)/2 Sqrt[2]], x/y}, {z, n Sqrt[5]}};
NCTeX[exp, DisplayPDF -> True, ImportPDF -> False]
uses DisplayPDF -> True
to ensure that the PDF viewer
is called and ImportPDF -> False
to prevent Mathematica
from displaying the formula inline. In other words, it displays the
formula in the PDF viewer without trying to import the PDF into
Mathematica. The default values for these options when using the
Mathematica notebook interface are:
DisplayPDF
(False
)ImportPDF
(True
)When NCTeX
is invoked using the command line interpreter
version of Mathematica the defaults are:
DisplayPDF
(False
)ImportPDF
(True
)Other useful options and their default options are:
Verbose
(False
),BreakEquations
(True
)TeXProcessor
(NCTeXForm
)Set BreakEquations -> True
to use the LaTeX package
beqn
to produce nice displays of long equations. Try the
following example:
expr = Series[Exp[x], {x, 0, 20}]
NCTeX[expr]
Use TexProcessor
to select your own TeX
converter. If NCAlgebra
is loaded then
NCTeXForm
is the default. Otherwise Mathematica’s
TeXForm
is used.
If Verbose -> True
you can see a detailed display of
what is going on behing the scenes. This is very useful for debugging.
For example, try:
expr = BesselJ[2, x]
NCTeX[exp, Verbose -> True]
to produce an output similar to the following one:
* NCTeX - LaTeX processor for NCAlgebra - Version 0.1
> Creating temporary file '/tmp/mNCTeX.tex'...
> Processing '/tmp/mNCTeX.tex'...
> Running 'latex -output-directory=/tmp/ /tmp/mNCTeX 1> "/tmp/mNCRun.out" 2> "/tmp/mNCRun.err"'...
> Running 'dvips -o /tmp/mNCTeX.ps -E /tmp/mNCTeX 1> "/tmp/mNCRun.out" 2> "/tmp/mNCRun.err"'...
> Running 'epstopdf /tmp/mNCTeX.ps 1> "/tmp/mNCRun.out" 2> "/tmp/mNCRun.err"'...
> Importing pdf file '/tmp/mNCTeX.pdf'...
Locate the files with extension .err as indicated by the verbose run of NCTeX to diagnose errors.
The remaining options:
PDFViewer
("open"
),LaTeXCommand
("latex"
)PDFLaTeXCommand
(Null
)DVIPSCommand
("dvips"
)PS2PDFCommand
("epstopdf"
)let you specify the names and, when appropriate, the path, of the
corresponding programs to be used by NCTeX
. Alternatively,
you can also directly implement custom versions of
NCRunDVIPS
NCRunLaTeX
NCRunPDFLaTeX
NCRunPDFViewer
NCRunPS2PDF
Those commands are invoked using NCRun
. Look at the
documentation for the package NCRun for more
details.
NCTeXForm
is a
replacement for Mathematica’s TeXForm
which adds
definitions allowing it to handle noncommutative expressions. It works
just as TeXForm
. NCTeXForm
is automatically
loaded with NCAlgebra
and is the default processor for
NCTeX
.
Here is an example:
SetNonCommutative[a, b, c, x, y];
exp = a ** x ** tp[b] - inv[c ** inv[a + b ** c] ** tp[y] + d]
NCTeXForm[exp]
produces
a.x.{b}^T-{\left(d+c.{\left(a+b.c\right)}^{-1}.{y}^T\right)}^{-1}
Note that the LaTeX output contains special code so that the
expression looks neat on the screen. You can see the result using
NCTeX
to convert the expression to PDF. Try
SetOptions[NCTeX, TeXProcessor -> NCTeXForm];
NCTeX[exp]
to produce
\(a.x.{b}^T-{\left(d+c.{\left (a+b.c\right)}^{-1}.{y}^T\right )}^{-1}\)
NCTeX
represents noncommutative products with a dot
(.
) in order to distinguish it from its commutative cousin.
We can see the difference in an expression that has both commutative and
noncommutative products:
exp = 2 a ** b - 3 c ** d
NCTeX[exp]
produces
\(2 \left(a.b\right) - 3 (c.d)\)
NCTeXForm handles lists and matrices as well. Here is a list:
exp = {x, tp[x], x + y, x + tp[y], x + inv[y], x ** x}
NCTeX[exp]
and its output:
\(\{ x, {x}^T, x+y, x+{y}^T, x+{y}^{-1}, x.x \}\)
and here is a matrix example:
exp = {{x, y}, {y, z}}
NCTeX[exp]
and its output:
\(\begin{bmatrix} x & y \\ y & z \end{bmatrix}\)
Here are some more examples:
exp = {{1 + Sin[x + (y - z)/2 Sqrt[2]], x/y}, {z, n Sqrt[5]}}
NCTeX[exp]
produces
\(\begin{bmatrix} 1+\operatorname{sin}{\left (x+\frac{1}{\sqrt{2}} \left (y-z\right )\right )} & x {y}^{-1} \\ z & \sqrt{5} n \end{bmatrix}\)
exp = {inv[x + y], inv[x + inv[y]]}
NCTeX[exp]
produces:
\(\{ {\left (x+y\right )}^{-1}, {\left (x+{y}^{-1}\right )}^{-1} \}\)
exp = {Sin[x], x y, Sin[x] y, Sin[x + y], Cos[gamma],
Sin[alpha] tp[x] ** (y - tp[y]), (x + tp[x]) (y ** z), -tp[y], 1/2,
Sqrt[2] x ** y}
NCTeX[exp]
produces:
\(\{ \operatorname{sin}{x}, x y, y \operatorname{sin}{x}, \operatorname{sin}{\left (x+y\right )}, \operatorname{cos}{\gamma}, \left({x}^T.\left (y-{y}^T\right )\right ) \operatorname{sin}{\alpha}, y z \left (x+{x}^T\right ), -{y}^T, \frac{1}{2}, \sqrt{2} \left(x.y\right ) \}\)
exp = inv[x + tp[inv[y]]]
NCTeX[exp]
produces:
\({\left (x+{{y}^T}^{-1}\right )}^{-1}\)
NCTeXForm
does not know as many functions as
TeXForm
. In some cases TeXForm
will produce
better results. Compare:
exp = BesselJ[2, x]
NCTeX[exp, TeXProcessor -> NCTeXForm]
output:
\(\operatorname{BesselJ}\left (2, x\right )\)
with
NCTeX[exp, TeXProcessor -> TeXForm]
output:
\(J_2(x)\)
It should be easy to customize NCTeXForm
though. Just
overload NCTeXForm
. In this example:
NCTeXForm[BesselJ[x_, y_]] := Format[BesselJ[x, y], TeXForm]
makes
NCTeX[exp, TeXProcessor -> NCTeXForm]
produce
\(J_2(x)\)
The following chapters and sections describes packages inside
NCAlgebra
. Detailed instructions for manual installation
are also provided in the section Manual
Installation.
Packages are automatically loaded unless otherwise noted.
Manual installation is no longer necessary nor advisable. As an alternative to manual installation, consider installing NCAlgebra via our paclet distribution, as discussed in section Installing NCAlgebra.
For those who do not use paclets, you can download NCAlgebra in one of the following ways.
git clone
You can clone the repository using git:
git clone https://github.com/NCAlgebra/NC
This will create a directory NC
which contains all files
neeeded to run NCAlgebra in any platform.
Cloning allows you to easily upgrade and switch between the various available releases. If you want to try the latest experimental version switch to branch devel using:
git checkout devel
If you’re happy with the latest stable release you do not need to do anything.
After you downloaded a zip file from github use your favorite zip
utility to unpack the file NC-master.zip
or
NC-devel.zip
on your favorite location.
IMPORTANT: Rename the top directory
NC
!
Releases are stable snapshots that you can find at
https://github.com/NCAlgebra/NC/releases
IMPORTANT: Rename the top directory
NC
!
Earlier releases can be downloaded from:
www.math.ucsd.edu/~ncalg
Releases in github are also tagged so you can easily switch from version to version using git.
All that is needed for NCAlgebra to run is that its top directory,
the NC
directory, be on Mathematica’s search path.
If you are on a unix flavored machine (Solaris, Linux, Mac OSX) then
unpacking or cloning in your home directory (~
) is all you
need to do.
Otherwise, you may need to add the installation directory to Mathematica’s search path.
If you are experienced with Mathematica:
Edit the main Mathematica init.m
file (not the
one inside the NC
directory) to add the name of the
directory which contains the NC
folder to the Mathematica
variable $Path
, as in:
AppendTo[$Path,"**YOUR_INSTALLATION_DIRECTORY**"];
You can locate your user init.m
file by typing:
FileNameJoin[{$UserBaseDirectory, "Kernel", "init.m"}]
in Mathematica.
NC is a meta package that enables the functionality of the NCAlgebra suite of non-commutative algebra packages for Mathematica.
If you installed the paclet version of NCAlgebra it is not necessary to load the context
NC
before loading otherNCAlgebra
packages.Loading the context
NC
in the paclet version is however still supported for backward compatibility. It does nothing more than posting the message:NC::Directory: You are using a paclet version of NCAlgebra.
The package can be loaded using Get
, as in
<< NC`
or Needs
, as in
Needs["NC`"]
Once NC
is loaded you will see a message like
NC::Directory: You are using the version of NCAlgebra which is found in: "/your_home_directory/NC".
or
NC::Directory: You are using a paclet version of NCAlgebra.
if you installed from our paclet distribution.
You can then proceed to load any other package from the NCAlgebra suite.
For example you can load the package NCAlgebra
using
<< NCAlgebra`
See section NCAlgebra and NCOptions for more options and details
available while loading NCAlgebra
.
The NC::Directory
message can be suppressed by using
standard Mathematica message control functions. For example,
Off[NC`NC::Directory]
<< NC`
or
Quiet[<< NC`, NC`NC::Directory]
will load NC
quietly. Note that you have to refer to the
message by its fully qualified name NC`NC::Directory
because the context NC
is only available after loading the
package.
NCAlgebra is the main package of the NCAlgebra suite of non-commutative algebra packages for Mathematica.
The package can be loaded using Get
, as in
<< NCAlgebra`
or Needs
, as in
Needs["NCAlgebra`"]
If the option SmallCapSymbolsNonCommutative
is
True
then NCAlgebra
will set all global single
letter small cap symbols as noncommutative. If that is not desired,
simply set SmallCapSymbolsNonCommutative
to
False
before loading NCAlgebra
, as in
<< NCOptions`
SetOptions[NCOptions, SmallCapSymbolsNonCommutative -> False];
<< NCAlgebra`
See NCOptions for details.
When loading NCAlgebra
, a message will be issued warning
users whether any letters have been set as noncommutative upon loading.
Those messages are documented here.
Users can use Mathematica’s Quiet
and Off
if
they do not want these messages to display. For example,
Off[NCAlgebra`NCAlgebra::SmallCapSymbolsNonCommutative]
<< NCAlgebra`
or
<< NCOptions`
SetOptions[NCOptions, SmallCapSymbolsNonCommutative -> False];
Off[NCAlgebra`NCAlgebra::NoSymbolsNonCommutative]
<< NCAlgebra`
will load NCAlgebra
without issuing a symbol assignment
message.
Upon loading NCAlgebra
for the first time, a large
banner will be shown. If you do not want this banner to be displayed set
the option ShowBanner
to
False
before loading, as in
<< NCOptions`
SetOptions[NCOptions, ShowBanner -> False];
<< NCAlgebra`
For example, the following commands will perform a completly quiet
loading of NC
and NCAlgebra
:
Quiet[<< NC`, NC`NC::Directory];
<< NCOptions`
SetOptions[NCOptions, ShowBanner -> False];
Quiet[<< NCAlgebra`, NCAlgebra`NCAlgebra::SmallCapSymbolsNonCommutative];
One of the following messages will be displayed after loading.
NCAlgebra::SmallCapSymbolsNonCommutative
, if small cap
single letter symbols have been set as noncomutative;NCAlgebra::NoSymbolsNonCommutative
, if no symbols have
been set as noncomutative by NCAlgebra
.The following options
can be set using
SetOptions
before loading other packages:
SmallCapSymbolsNonCommutative
(True
): If
True
, loading NCAlgebra
will set all global
single letter small cap symbols as noncommutative;ShowBanner
(True
): If True
, a
banner, when available, will be shown during the first loading of a
package.UseNotation
(False
): If True
use Mathematica’s package Notation
when setting pretty
output in NCSetOutput
.For example,
<< NC`
<< NCOptions`
SetOptions[NCOptions, ShowBanner -> False];
<< NCAlgebra`
suppress the NCAlgebra banner that is printed the first time NCAlgebra is loaded.
NonCommutativeMultiply
is the main package that provides
noncommutative functionality to Mathematica’s native
NonCommutativeMultiply
bound to the operator
**
.
Members are:
Aliases are:
aj[expr]
is the adjoint of expression expr
.
It is a conjugate linear involution.
co[expr]
is the conjugate of expression
expr
. It is a linear involution.
See also: aj.
Id
is noncommutative multiplicative identity. Actually
Id is now set equal 1
.
inv[expr]
is the 2-sided inverse of expression
expr
.
If Options[inv, Distrubute]
is False
(the
default) then
inv[a**b]
returns inv[a**a]
. Conversely, if
Options[inv, Distrubute]
is True
then it
returns inv[b]**inv[a]
.
rt[expr]
is the root of expression
expr
.
tp[expr]
is the tranpose of expression
expr
. It is a linear involution.
CommutativeQ[expr]
is True
if expression
expr
is commutative (the default), and False
if expr
is noncommutative.
See also: SetCommutative, SetNonCommutative.
NonCommutativeQ[expr]
is equal to
Not[CommutativeQ[expr]]
.
See also: CommutativeQ.
SetCommutative[a,b,c,...]
sets all the
Symbols
a
, b
, c
, …
to be commutative.
See also: SetNonCommutative, CommutativeQ, NonCommutativeQ.
SetCommutativeHold[a,b,c,...]
sets all the
Symbols
a
, b
, c
, …
to be commutative.
SetCommutativeHold
has attribute HoldAll
and can be used to set Symbols which have already been assigned a
value.
See also: SetNonCommutativeHold, SetCommutative, SetNonCommutative, CommutativeQ, NonCommutativeQ.
SetNonCommutative[a,b,c,...]
sets all the
Symbols
a
, b
, c
, …
to be noncommutative.
See also: SetCommutative, CommutativeQ, NonCommutativeQ.
SetNonCommutativeHold[a,b,c,...]
sets all the
Symbols
a
, b
, c
, …
to be noncommutative.
SetNonCommutativeHold
has attribute HoldAll
and can be used to set Symbols which have already been assigned a
value.
See also: SetCommutativeHold, SetCommutative, CommutativeQ, NonCommutativeQ.
SetCommutativeFunction[f]
sets expressions with
Head
f
, i.e. functions, to be commutative.
By default, expressions in which the Head
or any of its
arguments is noncommutative will be considered noncommutative. For
example,
SetCommutative[trace];
a ** b ** trace[a ** b]
evaluates to a ** b ** trace[a ** b]
while
SetCommutativeFunction[trace];
a ** b ** trace[a ** b]
evaluates to trace[a ** b] * a ** b
.
See also: SetCommutative, SetNonCommutative, CommutativeQ, NonCommutativeQ, tr.
SetNonCommutativeFunction[f]
sets expressions with
Head
f
, i.e. functions, to be non commutative.
This is only necessary if it has been previously set commutative by SetCommutativeFunction.
See also: SetCommutativeFunction, SetCommutative, SetNonCommutative, CommutativeQ, NonCommutativeQ, tr.
SNC
is an alias for SetNonCommutative
.
See also: SetNonCommutative.
SetCommutingOperators[a,b]
will define a rule that
substitute any noncommutative product b ** a
by
a ** b
, effectively making the pair a
and
b
commutative. If you want to create a rule to replace
a ** b
by b ** a
use
SetCommutingOperators[b,a]
instead.
See also: UnsetCommutingOperators, CommutingOperatorsQ
UnsetCommutingOperators[a,b]
remove any rules previously
created by SetCommutingOperators[a,b]
or
SetCommutingOperators[b,a]
.
See also: SetCommutingOperators, CommutingOperatorsQ
CommutingOperatorsQ[a,b]
returns True
if
a
and b
are commuting operators.
See also: SetCommutingOperators, UnsetCommutingOperators
NCNonCommutativeSymbolOrSubscriptQ[expr]
returns
True if expr
is an noncommutative symbol or a
noncommutative symbol subscript.
See also: NCNonCommutativeSymbolOrSubscriptExtendedQ, NCSymbolOrSubscriptQ, NCSymbolOrSubscriptExtendedQ, NCPowerQ.
NCNonCommutativeSymbolOrSubscriptExtendedQ[expr]
returns
True if expr
is an noncommutative symbol, a
noncommutative symbol subscript, or the transpose (tp
) or
adjoint (aj
) of a noncommutative symbol or noncommutative
symbol subscript.
See also: NCNonCommutativeSymbolOrSubscriptQ, NCSymbolOrSubscriptQ, NCSymbolOrSubscriptExtendedQ, NCPowerQ.
NCPowerQ[expr]
returns True if
expr
is an noncommutative symbol or symbol subscript or a
positive power of a noncommutative symbol or symbol subscript.
See also: NCNonCommutativeSymbolOrSubscriptQ, NCSymbolOrSubscriptQ.
Commutative[symbol]
is commutative even if
symbol
is noncommutative.
See also: CommuteEverything, CommutativeQ, SetCommutative, SetNonCommutative.
CommuteEverything[expr]
is an alias for BeginCommuteEverything.
See also: BeginCommuteEverything, Commutative.
BeginCommuteEverything[expr]
sets all symbols appearing
in expr
as commutative so that the resulting expression
contains only commutative products or inverses. It issues messages
warning about which symbols have been affected.
EndCommuteEverything[]
restores the symbols
noncommutative behaviour.
BeginCommuteEverything
answers the question what
does it sound like?
See also: EndCommuteEverything, Commutative.
EndCommuteEverything[expr]
restores noncommutative
behaviour to symbols affected by
BeginCommuteEverything
.
See also: BeginCommuteEverything, Commutative.
ExpandNonCommutativeMultiply[expr]
expands out
**
s in expr
.
For example
ExpandNonCommutativeMultiply[a**(b+c)]
returns
a**b + a**c.
NCExpand
is an alias for
ExpandNonCommutativeMultiply
.
See also: ExpandNonCommutativeMultiply, NCE.
NCE
is an alias for
ExpandNonCommutativeMultiply
.
See also: ExpandNonCommutativeMultiply, NCExpand.
NCExpandExponents[expr]
expands out powers of the
monomials appearing in expr
.
For example
NCExpandExponents[a**(b**c)^2**(c+d)]
returns
a**b**c**b**c**(c+d).
NCExpandExponents
only expands powers of monomials.
Powers of symbols or other expressions are not expanded using
NCExpandExponents
.
See also: NCToList ExpandNonCommutativeMultiply, NCExpand, NCE.
NCToList[expr]
produces a list with the symbols
appearing in monomial expr
. If expr
is not a
monomial it remains unevaluated. Powers of symbols are expanded before
the list is produced.
For example
NCToList[a**b**a^2]
returns
{a,b,a,a}
See also: NCExpandExponents, ExpandNonCommutativeMultiply, NCExpand, NCE.
NCTr
provides the commutative operator tr
that behaves as the standard mathematical trace operator.
Members are:
tr[expr]
is an linear operator with the following
properties:
tr
automatically distributes over sums;expr
is a noncommutative product, then product is
sorted; for example tr[b ** a]
evaluates into
tr[a ** b]
;tr[aj[expr]]
gets normalized as
Conjugate[tr[expr]]
.See also: SortCyclicPermutation, SortedCyclicPermutationQ.
SortedCyclicPermutation[list]
returns a cyclic
permutation of list sorted in ascending order.
See also: SortedCyclicPermutationQ.
SortCyclicPermutationQ[list]
returns True
if list
is a sorted cyclic permutation.
See also: SortCyclicPermutation.
Members are:
NCCollect[expr,vars]
collects terms of nc expression
expr
according to the elements of vars
and
attempts to combine them. It is weaker than NCStrongCollect in that only
same order terms are collected togther. It basically is
NCCompose[NCStrongCollect[NCDecompose]]]
.
If expr
is a rational nc expression then degree
correspond to the degree of the polynomial obtained using NCRationalToNCPolynomial.
NCCollect
also works with nc expressions instead of
Symbols in vars. In this case nc expressions are replaced by
new variables and NCCollect
is called using the resulting
expression and the newly created Symbols.
This command internally converts nc expressions into the special
NCPolynomial
format.
NCCollect[expr,vars,options]
uses options.
The following option is available:
ByTotalDegree
(False
): whether to collect
by total or partial degree.Notes:
While NCCollect[expr, vars]
always returns
mathematically correct expressions, it may not collect vars
from as many terms as one might think it should.
See also: NCStrongCollect, NCCollectSymmetric, NCCollectSelfAdjoint, NCStrongCollectSymmetric, NCStrongCollectSelfAdjoint, NCRationalToNCPolynomial.
NCCollectExponents[expr]
collects exponents in
noncommutative monomials.
For example
NCCollectExponents[a**b**a**b**c-a**b**a**b]
returns
(a**b)^2**c-(a**b)^2
See also: NCCollect, NCStrongCollect.
NCCollectSelfAdjoint[expr,vars]
allows one to collect
terms of nc expression expr
on the variables
vars
and their adjoints without writing out the
adjoints.
This command internally converts nc expressions into the special
NCPolynomial
format.
NCCollectSelfAdjoint[expr,vars,options]
uses
options.
The following option is available:
ByTotalDegree
(False
): whether to collect
by total or partial degree.See also: NCCollect, NCStrongCollect, NCCollectSymmetric, NCStrongCollectSymmetric, NCStrongCollectSelfAdjoint.
NCCollectSymmetric[expr,vars]
allows one to collect
terms of nc expression expr
on the variables
vars
and their transposes without writing out the
transposes.
This command internally converts nc expressions into the special
NCPolynomial
format.
NCCollectSymmetric[expr,vars,options]
uses options.
The following option is available:
ByTotalDegree
(False
): whether to collect
by total or partial degree.See also: NCCollect, NCStrongCollect, NCCollectSelfAdjoint, NCStrongCollectSymmetric, NCStrongCollectSelfAdjoint.
NCStrongCollect[expr,vars]
collects terms of expression
expr
according to the elements of vars
and
attempts to combine by association.
In the noncommutative case the Taylor expansion and so the collect
function is not uniquely specified. The function
NCStrongCollect
often collects too much and while correct
it may be stronger than you want.
For example, a symbol x
will factor out of terms where
it appears both linearly and quadratically thus mixing orders.
This command internally converts nc expressions into the special
NCPolynomial
format.
See also: NCCollect, NCCollectSymmetric, NCCollectSelfAdjoint, NCStrongCollectSymmetric, NCStrongCollectSelfAdjoint.
NCStrongCollectSymmetric[expr,vars]
allows one to
collect terms of nc expression expr
on the variables
vars
and their transposes without writing out the
transposes.
This command internally converts nc expressions into the special
NCPolynomial
format.
See also: NCCollect, NCStrongCollect, NCCollectSymmetric, NCCollectSelfAdjoint, NCStrongCollectSymmetric.
NCStrongCollectSymmetric[expr,vars]
allows one to
collect terms of nc expression expr
on the variables
vars
and their transposes without writing out the
transposes.
This command internally converts nc expressions into the special
NCPolynomial
format.
See also: NCCollect, NCStrongCollect, NCCollectSymmetric, NCCollectSelfAdjoint, NCStrongCollectSelfAdjoint.
NCCompose[dec]
will reassemble the terms in
dec
which were decomposed by NCDecompose
.
NCCompose[dec, degree]
will reassemble only the terms of
degree degree
.
The expression NCCompose[NCDecompose[p,vars]]
will
reproduce the polynomial p
.
The expression NCCompose[NCDecompose[p,vars], degree]
will reproduce only the terms of degree degree
.
This command internally converts nc expressions into the special
NCPolynomial
format.
See also: NCDecompose, NCPDecompose.
NCDecompose[p,vars]
gives an association of elements of
the nc polynomial p
in variables vars
in which
elements of the same order are collected together.
NCDecompose[p]
treats all nc letters in p
as variables.
This command internally converts nc expressions into the special
NCPolynomial
format.
Internally NCDecompose
uses NCPDecompose
.
See also: NCCompose, NCPDecompose.
NCTermsOfDegree[expr,vars,degrees]
returns an expression
such that each term has degree degrees
in variables
vars
.
For example,
NCTermsOfDegree[x**y**x - x**x**y + x**w + z**w, {x,y}, {2,1}]
returns x**y**x - x**x**y
,
NCTermsOfDegree[x**y**x - x**x**y + x**w + z**w, {x,y}, {1,0}]
returns x**w
,
NCTermsOfDegree[x**y**x - x**x**y + x**w + z**w, {x,y}, {0,0}]
returns z**w
, and
NCTermsOfDegree[x**y**x - x**x**y + x**w + z**w, {x,y}, {0,1}]
returns 0
.
This command internally converts nc expressions into the special
NCPolynomial
format.
See also: NCTermsOfTotalDegree, NCDecompose, NCPDecompose.
NCTermsOfTotalDegree[expr,vars,degree]
returns an
expression such that each term has total degree degree
in
variables vars
.
For example,
NCTermsOfTotalDegree[x**y**x - x**x**y + x**w + z**w, {x,y}, 3]
returns x**y**x - x**x**y
,
NCTermsOfTotalDegree[x**y**x - x**x**y + x**w + z**w, {x,y}, 1]
returns x**w
,
NCTermsOfTotalDegree[x**y**x - x**x**y + x**w + z**w, {x,y}, 0]
returns z**w
, and
NCTermsOfTotalDegree[x**y**x - x**x**y + x**w + z**w, {x,y}, 2]
returns 0
.
This command internally converts nc expressions into the special
NCPolynomial
format.
See also: NCTermsOfDegree, NCDecompose, NCPDecompose.
NCReplace is a package containing several functions
that are useful in making replacements in noncommutative expressions. It
offers replacements to Mathematica’s Replace
,
ReplaceAll
, ReplaceRepeated
, and
ReplaceList
functions.
Commands in this package replace the old
Substitute
andTransform
family of command which are been deprecated. The new commands are much more reliable and work faster than the old commands. From the beginning, substitution was always problematic and certain patterns would be missed. We reassure that the call expression that are returned are mathematically correct but some opportunities for substitution may have been missed.
Members are:
Aliases:
Options:
ApplyPowerRule
(True
): If
True
, NCReplacePowerRule
is automatically
applied to all rules in NCReplace
,
NCReplaceAll
, NCReplaceRepeated
, and
NCReplaceList
.NCReplace[expr,rules]
applies a rule or list of rules
rules
in an attempt to transform the entire nc expression
expr
.
NCReplace[expr,rules,levelspec]
applies
rules
to parts of expr
specified by
levelspec
.
NCReplace[expr,rules,options]
and
NCReplace[expr,rules,levelspec]
use
options
.
See also: NCReplaceAll, NCReplaceList, NCReplaceRepeated.
NCReplaceAll[expr,rules]
applies a rule or list of rules
rules
in an attempt to transform each part of the nc
expression expr
.
NCReplaceAll[expr,rules,options]
uses
options
.
See also: NCReplace, NCReplaceList, NCReplaceRepeated.
NCReplace[expr,rules]
attempts to transform the entire
nc expression expr
by applying a rule or list of rules
rules
in all possible ways, and returns a list of the
results obtained.
ReplaceList[expr,rules,n]
gives a list of at most
n
results.
NCReplaceList[expr,rules,options]
and
NCReplaceList[expr,rules,n,options]
use
options
.
See also: NCReplace, NCReplaceAll, NCReplaceRepeated.
NCReplaceRepeated[expr,rules]
repeatedly performs
replacements using rule or list of rules rules
until
expr
no longer changes.
NCReplaceRepeated[expr,rules,options]
uses
options
.
See also: NCReplace, NCReplaceAll, NCReplaceList, NCExpandReplaceRepeated.
NCExpandReplaceRepeated[expr,rules]
repeatedly applies
NCReplaceRepeated
and NCExpand
until the
results does not change.
NCExpandReplaceRepeated[expr,rules,options]
uses
options
.
See also: NCReplaceRepeated, NCExpandReplaceRepeatedSymmetric, NCExpandReplaceRepeatedSelfAdjoint.
NCR
is an alias for NCReplace
.
See also: NCReplace.
NCRA
is an alias for NCReplaceAll
.
See also: NCReplaceAll.
NCRR
is an alias for NCReplaceRepeated
.
See also: NCReplaceRepeated.
NCRL
is an alias for NCReplaceList
.
See also: NCReplaceList.
NCMakeRuleSymmetric[rules]
add rules to transform the
transpose of the left-hand side of rules
into the transpose
of the right-hand side of rules
.
See also: NCMakeRuleSelfAdjoint, NCReplace, NCReplaceAll, NCReplaceList, NCReplaceRepeated.
NCMakeRuleSelfAdjoint[rules]
add rules to transform the
adjoint of the left-hand side of rules
into the adjoint of
the right-hand side of rules
.
See also: NCMakeRuleSymmetric, NCReplace, NCReplaceAll, NCReplaceList, NCReplaceRepeated.
NCReplaceSymmetric[expr, rules]
applies
NCMakeRuleSymmetric
to rules
before calling
NCReplace
.
NCReplaceSymmetric[expr,rules,options]
uses
options
.
See also: NCReplace, NCMakeRuleSymmetric.
NCReplaceAllSymmetric[expr, rules]
applies
NCMakeRuleSymmetric
to rules
before calling
NCReplaceAll
.
NCReplaceAllSymmetric[expr,rules,options]
uses
options
.
See also: NCReplaceAll, NCMakeRuleSymmetric.
NCReplaceRepeatedSymmetric[expr, rules]
applies
NCMakeRuleSymmetric
to rules
before calling
NCReplaceRepeated
.
NCReplaceRepeatedSymmetric[expr,rules,options]
uses
options
.
See also: NCReplaceRepeated, NCMakeRuleSymmetric.
NCExpandReplaceRepeatedSymmetric[expr,rules]
repeatedly
applies NCReplaceRepeatedSymmetric
and
NCExpand
until the results does not change.
NCExpandReplaceRepeatedSymmetric[expr,rules,options]
uses options
.
See also: NCReplaceRepeatedSymmetric, NCExpandReplaceRepeated, NCExpandReplaceRepeatedSelfAdjoint.
NCReplaceListSymmetric[expr, rules]
applies
NCMakeRuleSymmetric
to rules
before calling
NCReplaceList
.
NCReplaceListSymmetric[expr,rules,options]
uses
options
.
See also: NCReplaceList, NCMakeRuleSymmetric.
NCRSym
is an alias for
NCReplaceSymmetric
.
See also: NCReplaceSymmetric.
NCRASym
is an alias for
NCReplaceAllSymmetric
.
See also: NCReplaceAllSymmetric.
NCRRSym
is an alias for
NCReplaceRepeatedSymmetric
.
See also: NCReplaceRepeatedSymmetric.
NCRLSym
is an alias for
NCReplaceListSymmetric
.
See also: NCReplaceListSymmetric.
NCReplaceSelfAdjoint[expr, rules]
applies
NCMakeRuleSelfAdjoint
to rules
before calling
NCReplace
.
NCReplaceSelfAdjoint[expr,rules,options]
uses
options
.
See also: NCReplace, NCMakeRuleSelfAdjoint.
NCReplaceAllSelfAdjoint[expr, rules]
applies
NCMakeRuleSelfAdjoint
to rules
before calling
NCReplaceAll
.
NCReplaceAllSelfAdjoint[expr,rules,options]
uses
options
.
See also: NCReplaceAll, NCMakeRuleSelfAdjoint.
NCReplaceRepeatedSelfAdjoint[expr, rules]
applies
NCMakeRuleSelfAdjoint
to rules
before calling
NCReplaceRepeated
.
NCReplaceRepeatedSelfAdjoint[expr,rules,options]
uses
options
.
See also: NCReplaceRepeated, NCMakeRuleSelfAdjoint.