SciSharp/Numpy.NET

got err when i use 'np.linalg.qr'

bigpo opened this issue · 2 comments

bigpo commented

hi, guys, Could anybody show me how to use the func np.linalg.qr ?
when i use it like this

var a = np.random.randn(3, 3);
var tmp = np.linalg.qr(a);

i got errors:

Unhandled exception. Python.Runtime.PythonException: IndexError : tuple index out of range
at Python.Runtime.PyObject.GetItem(PyObject key)
at Python.Runtime.PyObject.GetItem(Int32 index)
at Python.Runtime.PyObject.get_Item(Int32 index)
at System.Dynamic.UpdateDelegates.UpdateAndExecute2[T0,T1,TRet](CallSite site, T0 arg0, T1 arg1)
at Numpy.np.linalg.qr(NDarray a, String mode)
at ConsoleApp2.Program.Main(String[] args) in C:\Users\xx\source\repos\ConsoleApp2\Program.cs:line 11 

Thanks!

henon commented

I fixed qr. It returns different things depending on the parameters it gets. C# can not do the same as python here. So you have to take either the first, the first and second or all three results of the function depending on which mode you use. In the cases where python returns only one argument the other two entries of the tuple are null. you can see how to use it in this unit test:

        [TestMethod]
        public void qrTest()
        {
            // >>> a = np.random.randn(9, 6)
            // >>> q, r = np.linalg.qr(a)
            // >>> np.allclose(a, np.dot(q, r))  # a does equal qr
            // True
            // >>> r2 = np.linalg.qr(a, mode='r')
            // >>> r3 = np.linalg.qr(a, mode='economic')
            // >>> np.allclose(r, r2)  # mode='r' returns the same r as mode='full'
            // True
            // >>> # But only triu parts are guaranteed equal when mode='economic'
            // >>> np.allclose(r, np.triu(r3[:6,:6], k=0))
            // True
            // 

            {
                var a = np.random.randn(9, 6);
                var (q, r, _) = np.linalg.qr(a);
                //Console.WriteLine("r: " + r.repr);
                var given = np.allclose(a, np.dot(q, r)); // a does equal qr;
                Assert.AreEqual(true, given);
                (var r2, _, _) = np.linalg.qr(a, mode: "r");
                //Console.WriteLine("r2: " + r2.repr);
                (var r3, _, _) = np.linalg.qr(a, mode: "economic");
                given = np.allclose(r, r2); // mode='r' returns the same r as mode='full';
                Assert.AreEqual(true, given);
            }
#if TODO
            // But only triu parts are guaranteed equal when mode='economic';
            given=  np.allclose(r, np.triu(r3{:6,:6}, k=0));
             expected=
                "True";
            Assert.AreEqual(expected, given.repr);
#endif 

            // Example illustrating a common use of qr: solving of least squares
            // problems

            // What are the least-squares-best m and y0 in y = y0 + mx for
            // the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
            // and you’ll see that it should be y0 = 0, m = 1.)  The answer is provided
            // by solving the over-determined matrix equation Ax = b, where:

            // A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
            // x = array([[y0], [m]])
            // b = array([[1], [0], [2], [1]])
            // 

            // If A = qr such that q is orthonormal (which is always possible via
            // Gram-Schmidt), then x = inv(r) * (q.T) * b.  (In numpy practice,
            // however, we simply use lstsq.)

            // >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
            // >>> A
            // array([[0, 1],
            //        [1, 1],
            //        [1, 1],
            //        [2, 1]])
            // >>> b = np.array([1, 0, 2, 1])
            // >>> q, r = LA.qr(A)
            // >>> p = np.dot(q.T, b)
            // >>> np.dot(LA.inv(r), p)
            // array([  1.1e-16,   1.0e+00]) // <--- note henon: this result seems to be wrong in the documentation.
            // 
            {
                var A = np.array(new[,] { { 0, 1 }, { 1, 1 }, { 1, 1 }, { 2, 1 } });
                NDarray given = A;
                var expected =
                    "array([[0, 1],\n" +
                    "       [1, 1],\n" +
                    "       [1, 1],\n" +
                    "       [2, 1]])";
                Assert.AreEqual(expected, given.repr);
                var b = np.array(new[]{ 1, 0, 2, 1 });
                (var q, var r, _) = np.linalg.qr(A);
                //Console.WriteLine("q:" + q.repr);
                //Console.WriteLine("r:"+r.repr);
                var p = np.dot(q.T, b);
                //Console.WriteLine("p:" + p.repr);
                var r_inv = np.linalg.inv(r);
                //Console.WriteLine("r_inv:" + r_inv.repr);
                given = np.dot(r_inv, p);
                expected =
                    "array([6.66133815e-16, 1.00000000e+00])";
                Assert.AreEqual(expected, given.repr);
            }
        }

The test runs through now, so it should work now. Publishing the fix to nuget in a few minutes

bigpo commented

It works, Thank you very much!