MichaĆ«l Baudin [Fri, 8 Mar 2013 16:21:37 +0000 (17:21 +0100)]
Change-Id: I2873051ae08f0b591f3dcbc7a50c9d80d03123e0

index 57b8719..1a5baff 100644 (file)
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
-        <synopsis>x=A\b</synopsis>
+        <synopsis>x=A\B</synopsis>
</refsynopsisdiv>
<refsection>
<title>Description</title>
<para>
-            Backslash denotes left matrix division.
-            <literal>x=A\b</literal> is a solution to <literal>A*x=b</literal>.
+            Backslash is the left matrix division:
+            <literal>x=A\B</literal> is a solution to <literal>A*x=B</literal>.
</para>
<para>
-            If <literal>A</literal> is square and non-singular <literal>x=A\b</literal> (uniquely defined) is equivalent to <literal>x=inv(A)*b</literal> (but the computations are much cheaper).
+            If <literal>A</literal> is square and non-singular <literal>x=A\B</literal> is
+            equivalent to <literal>x=inv(A)*B</literal> in exact arithmetic,
+            but the computations are more accurate and cheaper in floating point arithmetic.
+            Hence, to compute the solution of the linear system of equations <literal>A*x=B</literal>,
+            the backslash operator should be used, and the <literal>inv</literal> function
+            should be avoided.
+        </para>
+        <para>
+            In the case where A is square, the solution <literal>x</literal> can be computed
+            either from LU factorization or from a linear least squares solver.
+            If the condition number of A is smaller than 1/(10*%eps) (i.e. if A is well conditionned),
+            the LU factorization with row pivoting is used.
+            If not (i.e. if A is poorly conditionned), then X is the minimum-norm solution which
+            minimizes <literal>||A*X-B||</literal> using a complete
+            orthogonal factorization of A (i.e. X is the solution of a linear least squares problem).
</para>
<para>
If <literal>A</literal> is not square, <literal>x</literal> is a least square solution,
-            i.e. <literal>norm(A*x-b)</literal> is minimal (Euclidean norm). If <literal>A</literal> is full
-            column rank, the least square solution, <literal>x=A\b</literal>, is uniquely
-            defined (there is a unique <literal>x</literal> which minimizes <literal>norm(A*x-b)</literal>).
+            i.e. <literal>norm(A*x-B)</literal> is minimal (Euclidean norm). If <literal>A</literal> is full
+            column rank, the least square solution, <literal>x=A\B</literal>, is uniquely
+            defined (there is a unique <literal>x</literal> which minimizes <literal>norm(A*x-B)</literal>).
If <literal>A</literal> is not full column rank, then the least square
-            solution is not unique, and <literal>x=A\b</literal>, in general, is not the solution
-            with minimum norm (the minimum norm solution is <literal>x=pinv(A)*b</literal>).
+            solution is not unique, and <literal>x=A\B</literal>, in general, is not the solution
+            with minimum norm (the minimum norm solution is <literal>x=pinv(A)*B</literal>).
</para>
<para>
<literal>A.\B</literal>  is the matrix with <literal>(i,j)</literal> entry <literal>A(i,j)\B(i,j)</literal>.
-            If <literal>A</literal> (or <literal>B</literal>) is a scalar <literal>A.\B</literal> is equivalent to
+            If <literal>A</literal> (or <literal>B</literal>) is a scalar <literal>A.\B</literal> is equivalent to
<literal>A*ones(B).\B</literal> (or <literal>A.\(B*ones(A))</literal>.
</para>
<para>
<literal>A\.B</literal>  is an operator with no predefined meaning. It may be used
+            the same precedence as <literal>*</literal> or <literal>/</literal>.
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
-A=rand(3,2);b=[1;1;1]; x=A\b; y=pinv(A)*b;  x-y
-A=rand(2,3);b=[1;1]; x=A\b; y=pinv(A)*b; x-y, A*x-b, A*y-b
+A=[
+   9.   -36.    30.
+  -36.   192.  -180.
+   30.  -180.   180.
+];
+b=[
+   3.
+  -24.
+   30.
+];
+x=A\b
+A*x-b // close to zero
+
+A=rand(3,2);
+b=[1;1;1];
+x=A\b;
+y=pinv(A)*b;
+x-y
+A=rand(2,3);b=[1;1];
+x=A\b;
+y=pinv(A)*b;
+x-y, A*x-b, A*y-b

// if rank is deficient
-A=rand(3,1)*rand(1,2); b=[1;1;1]; x=A\b; y=pinv(A)*b; A*x-b, A*y-b
-A=rand(2,1)*rand(1,3); b=[1;1]; x=A\b; y=pinv(A)*b; A*x-b, A*y-b
+A=rand(3,1)*rand(1,2);
+b=[1;1;1];
+x=A\b;
+y=pinv(A)*b;
+A*x-b, A*y-b
+A=rand(2,1)*rand(1,3);
+b=[1;1];
+x=A\b;
+y=pinv(A)*b;
+A*x-b, A*y-b

// A benchmark of several linear solvers

+   "/modules/umfpack/examples/bcsstk24.rsa");

b = zeros(size(A,1),1);

tic();
res = umfpack(A,'\',b);
-mprintf('\ntime needed to solve the system with umfpack: %.3f\n',toc());
+mprintf('\nTime with umfpack: %.3f\n',toc());

tic();
res = linsolve(A,b);
-mprintf('\ntime needed to solve the system with linsolve: %.3f\n',toc());
+mprintf('\ntime with linsolve: %.3f\n',toc());

tic();
res = A\b;
-mprintf('\ntime needed to solve the system with the backslash operator: %.3f\n',toc());
+mprintf('\ntime with backslash: %.3f\n',toc());
]]></programlisting>
</refsection>
<refsection role="see also">
@@ -91,4 +135,13 @@ mprintf('\ntime needed to solve the system with the backslash operator: %.3f\n',
</member>
</simplelist>
</refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.4.1</revnumber>
+                <revremark>The threshold level for conditioning in blackslash increased.</revremark>
+            </revision>
+        </revhistory>
+    </refsection>
</refentry>
index 7cd3096..343a50e 100644 (file)
@@ -2,31 +2,40 @@
<refnamediv>
<refname>slash</refname>
-        <refpurpose>(/) right division and feed back (comments)</refpurpose>
+        <refpurpose>(/) right division and feed back</refpurpose>
</refnamediv>
<refsection>
<title>Description</title>
<para>
-            Right division.  <code>x = A / b</code> is the solution of <code>x * b = A</code>.
+            Right division: <code>x=A/B</code> is the solution of <code>x*B=A</code>.
</para>
<para>
-            <code>b/a = (a' \ b')'</code>.
+            The slash (right division) and backslash (left division) operators are linked by the following equation:
+            <code>B/A=(A'\B')'</code>.
</para>
<para>
-            <code>a ./ b</code> is the matrix with entries <literal>a(i,j)/ b(i,j)</literal>.
-            If <literal>b</literal> is scalar (1x1 matrix) this operation is the same
-            as <code> a ./ b * ones(a)</code>. (Same convention if <literal>a</literal> is a scalar).
+            In the case where A is square, the solution <literal>x</literal> can be computed
+            either from LU factorization or from a linear least squares solver.
+            If the condition number of A is smaller than 1/(10*%eps) (i.e. if A is well conditionned),
+            the LU factorization with row pivoting is used.
+            If not (i.e. if A is poorly conditionned), then X is the minimum-norm solution which
+            minimizes <literal>||A*X-B||</literal> using a complete
+            orthogonal factorization of A (i.e. X is the solution of a linear least squares problem).
+        </para>
+        <para>
+            <code>A./B</code> is the elementwise right division, i.e. the matrix with
+            entries <literal>A(i,j)/B(i,j)</literal>.
+            If <literal>B</literal> is scalar (1x1 matrix) this operation is the same
+            as <code>A./B*ones(A)</code>.
+            Same convention if <literal>A</literal> is a scalar.
</para>
<para>
<note>
-                Remark that <code>123./b</code> is interpreted as <code>(123.)/b</code>. In this
+                Remark that <code>123./B</code> is interpreted as <code>(123.)/B</code>. In this
cases dot is part of the  number not of the operator.
</note>
</para>
<para>
-            Backslash stands for left division.
-        </para>
-        <para>
System feed back.  <code>S = G/.K</code> evaluates
<code>S = G*(eye() + K*G)^(-1)</code> this operator avoid simplification problem.
</para>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
+a=[3.,-24.,30.];
+B=[
+   9.   -36.    30.
+  -36.   192.  -180.
+   30.  -180.   180.
+];
+x=a/B
+x*B-a // close to zero
+
a=4 / 2; // Should be 2
a=2 ./ [2,4]; //     1.    0.5
// Comments are good. They help to understand code
@@ -63,4 +81,14 @@ a=2 ./ [2,4]; //     1.    0.5
</member>
</simplelist>
</refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.4.1</revnumber>
+                <revremark>The threshold level for conditioning in slash increased.</revremark>
+            </revision>
+        </revhistory>
+    </refsection>
+
</refentry>