{"action":"created","issue":{"url":"https://api.github.com/repos/arrayfire/arrayfire/issues/2143","repository_url":"https://api.github.com/repos/arrayfire/arrayfire","labels_url":"https://api.github.com/repos/arrayfire/arrayfire/issues/2143/labels{/name}","comments_url":"https://api.github.com/repos/arrayfire/arrayfire/issues/2143/comments","events_url":"https://api.github.com/repos/arrayfire/arrayfire/issues/2143/events","html_url":"https://github.com/arrayfire/arrayfire/issues/2143","id":318535274,"node_id":"MDU6SXNzdWUzMTg1MzUyNzQ=","number":2143,"title":"Least squares","user":{"login":"FloopCZ","id":6018335,"node_id":"MDQ6VXNlcjYwMTgzMzU=","avatar_url":"https://avatars1.githubusercontent.com/u/6018335?v=4","gravatar_id":"","url":"https://api.github.com/users/FloopCZ","html_url":"https://github.com/FloopCZ","followers_url":"https://api.github.com/users/FloopCZ/followers","following_url":"https://api.github.com/users/FloopCZ/following{/other_user}","gists_url":"https://api.github.com/users/FloopCZ/gists{/gist_id}","starred_url":"https://api.github.com/users/FloopCZ/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/FloopCZ/subscriptions","organizations_url":"https://api.github.com/users/FloopCZ/orgs","repos_url":"https://api.github.com/users/FloopCZ/repos","events_url":"https://api.github.com/users/FloopCZ/events{/privacy}","received_events_url":"https://api.github.com/users/FloopCZ/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":5,"created_at":"2018-04-27T20:11:14Z","updated_at":"2018-07-06T20:10:46Z","closed_at":null,"author_association":"CONTRIBUTOR","body":"Hi everyone,\r\n\r\nI would like to know what is the \"official\" way of calculating least squares approximation of a system of linear equations, i.e., `Ax = b` for non-square `A`.\r\n\r\nUp to this moment, I believed that when `af::solve` gets an overdetermined system of equations (e.g., `af::solve(af::randu(4, 2), af::randu(4))`), it produces the least squares approximation of the coefficients (because in the documentation, it is in a section called [Linear solve and least squares](http://arrayfire.org/docs/group__lapack__solve__func__gen.htm), although in the function description, no such thing is mentioned). Unfortunately, I just performed some tests and it is not the case. My questions are:\r\n\r\n1. Is there a recommended way to solving an overdetermined system of linear equations?\r\n2. What do the numbers computed by `af::solve` mean? Shouldn't it rather throw an exception if this functionality is not supported?\r\n\r\nNotes:\r\n- `af::inverse` does not work for non-squared matrices.\r\n- Relates to #2074 and #2077.\r\n- Thanks @Myhroun for detecting the issue in the first place.\r\n\r\nThank you and have a nice weekend!"},"comment":{"url":"https://api.github.com/repos/arrayfire/arrayfire/issues/comments/403134587","html_url":"https://github.com/arrayfire/arrayfire/issues/2143#issuecomment-403134587","issue_url":"https://api.github.com/repos/arrayfire/arrayfire/issues/2143","id":403134587,"node_id":"MDEyOklzc3VlQ29tbWVudDQwMzEzNDU4Nw==","user":{"login":"FloopCZ","id":6018335,"node_id":"MDQ6VXNlcjYwMTgzMzU=","avatar_url":"https://avatars1.githubusercontent.com/u/6018335?v=4","gravatar_id":"","url":"https://api.github.com/users/FloopCZ","html_url":"https://github.com/FloopCZ","followers_url":"https://api.github.com/users/FloopCZ/followers","following_url":"https://api.github.com/users/FloopCZ/following{/other_user}","gists_url":"https://api.github.com/users/FloopCZ/gists{/gist_id}","starred_url":"https://api.github.com/users/FloopCZ/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/FloopCZ/subscriptions","organizations_url":"https://api.github.com/users/FloopCZ/orgs","repos_url":"https://api.github.com/users/FloopCZ/repos","events_url":"https://api.github.com/users/FloopCZ/events{/privacy}","received_events_url":"https://api.github.com/users/FloopCZ/received_events","type":"User","site_admin":false},"created_at":"2018-07-06T20:10:45Z","updated_at":"2018-07-06T20:10:45Z","author_association":"CONTRIBUTOR","body":"Here is a bit more elegant and vectorized code for least squares and pseudoinverse with some additional comments.\r\n```c++\r\n// Pseudoinverse matrix.\r\naf::array pinv(const af::array& a)\r\n{\r\n    af::array u, s, vt;\r\n    af::svd(u, s, vt, a);\r\n    // Having SVD decomposition, we need the inverse of the diagonal matrix S.\r\n    // That is, a diagonal matrix with 1/S on its main diagonal. On the positions\r\n    // where S is (close to) 0, this inverse matrix is set to 0 as well.\r\n    af::array s_pinv = af::constant(0, a.dims(0), a.dims(1), a.type());\r\n    af::seq diagonal(0, (s.elements() - 1) * (s_pinv.dims(0) + 1), s_pinv.dims(0) + 1);\r\n    // The S is considered to be 0 when it is smaller than the following threshold:\r\n    double threshold = 1e-6 * std::max(a.dims(0), a.dims(1)) * af::max<double>(s);\r\n    s_pinv(diagonal) = af::select(s < threshold, 0., 1. / s);\r\n    // The final pseudoinverse is V * S^-1 * U^T.\r\n    return af::matmul(af::matmul(vt.T(), s_pinv.T()), u.T());\r\n}\r\n\r\n// Get x, such that a*x=b minimizes the mean squared error.\r\naf::array lstsq(const af::array& a, const af::array& b)\r\n{\r\n    return af::matmul(pinv(a), b);\r\n}\r\n\r\n// Usage:\r\naf::array x = lstsq(a, b);\r\n```"}}
{"actor":{"display_login":"FloopCZ"}}